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FOREWORD 


This  report  describes  the  work  performed  by  Bell  Aerospace  Textron, 
a  Division  of  Textron,  Inc.  Buffalo,  New  York.  The  work  was  sponsored  by 
the  Flight  Dynamics  Laboratory,  Air  Force  Wright  Aeronautical  Laboratories, 
Wright-Patterson  Air  Force  Base,  Ohio,  under  Contract  F33615-80-C3214. 

The  work  was  initiated  under  Project  2307,  "Research  in  Flight 
Vehicle  Dynamics"  Task  2307N518,  "Basic  Research  in  Structures  and  Dynamics". 
The  work  was  administered  by  Dr.  N.  S.  Khot,  Project  Engineer  of  the 
Structures  and  Dynamics  Division  (FIBRA). 

The  contracted  work  was  performed  between  August  1980  and  December 

1982. 

The  work  was  performed  in  the  Structures  and  Vehicles  Systems 
Directorate,  Bell  Aerospace  Textron.  Mr.  James  R.  Batt  was  the  Program 
Manager /Technical  Director  of  the  study. 
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1.0  INTRODUCTION 


Modern  technology  of  structural  optimization  is  now  nearly  two 
decades  old.  During  that  period  an  intense  amount  of  effort  has  been 
expended  on  studying  the  very  many  facets  of  the  optimization  problem.  At 
a  conservative  estimate  a  bibliography  of  over  300  references  to  relevant 
work  in  this  area  could  be  compiled  and,  in  fact,  the  search  for  the  most 

* 

efficient  structure  to  perform  a  given  task  is  much  older  than  20  years. 

«  The  standard  references  to  Michell,  Ref.  1,  and  Clerk-Maxwell,  Ref.  2, 

can  be  made  to  illustrate  that  the  ideals  and  objectives  of  designing 
practical  efficient  structures  are  very  fertile  fields  which  have  been 
studied  in  the  extreme.  Yet,  in  comparison  with  the  erstwhile  companion 
technology  of  finite  element  analysis,  it  must  be  recognized,  regretfully, 
that  structural  optimization  methodology  has  had  little,  if  any,  real  im¬ 
pact  on  the  design  of  modern  structural  systems.  It  has  simply  been  not 
possible,  in  general,  to  translate  the  vast  research  into  acceptable  prac¬ 
tical  design  tools.  While  a  number  of  operational  optimization  programs 
have  been  developed  -  of  various  types  -  their  use  and  acceptance  has  been 
very  limited.  There  are  exceptions,  but  there  has  been  no  general  utiliza¬ 
tion  of  any  optimization  methods  such  as  has  occurred  in  the  universal  use 

r 

of  finite  element  technology. 

'  The  reasons  for  this  are  manifold.  Firstly,  optimization  represents 

a  degree  of  sophistication  above  pure  analysis  and  hence  is  more  costly. 
Secondly,  the  sophistication  of  mathematics  may  be  suspect  by  designers 
and  engineers  who  have  no  good  reference  for  judging  the  validity  of  the 
concepts.  Thirdly,  the  lack  of  a  unified  thrust  on  the  part  of  researchers 
and  developers  to  select  an  approach  to  the  problem  is  strongly  indicative 
of  confidence  in  the  value  and  validity  of  structural  optimization  and 
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fourthly,  the  nagging  concern  of  deaignera  about  the  "goodness"  of  an 
optimized  structure.  This  final  point  is  one  which  has  long  been  the  cen¬ 
tral  point  of  a  major  controversy.  The  removal  of  concerns  about  the  re¬ 
liability  of  optimized  structures  should  provide  a  major  step  forward  in 
the  widespread  acceptance  of  optimization  technology. 

In  this  context,  it  is  most  relevant  to  consider  the  appropriate 
state-of-the-art  in  structural  optimization  and  its  applicability  to  the 
present  study  reported  herein. 

In  spite  of  the  very  large  effort  of  work  expended  by  so  many  re¬ 
searchers,  many  of  the  most  basic  problems  are  still  largely  unsolved. 
Perhaps  one  of  the  most  significant  areas  of  progress  has  been  in  scoping 
our  domain  of  ignorance.  From  the  early  days,  when  it  was  first  proposed 
that  some  methods  of  operations  research  could  be  linked  with  iterative 
finite  element  analysis  to  the  present  day  studies  using  Lagrangian 
formulations,  there  has  been  a  much  greater  understanding  of  the  true 
nature  of  the  problems  to  be  solved.  The  problems  themselves  have  not 
necessarily  been  solved  -  they  have  just  been  identified  more  precisely. 

In  structural  optimization,  with  the  usual  limitations  of  fixed 
geometry,  etc.,  there  is  obviously  a  hierarchy  of  relevant  constraints  - 
number  stresses,  nodal  deflections,  local  and  overall  stability,  frequency 
response,  flutter,  etc.  The  precise  order  of  higher  order  members  in  the 
hierarchy  is  more  subjective  and  is  influenced  by  the  types  of  structures 
under  consideration.  The  first  two,  strength  and  stiffness  are,  generally, 
prime  requirements  in  virtually  every  structure. 

Similarly,  there  are,  theoretically,  a  number  of  potential  candidate 
figures  of  merit  by  which  the  optimality  of  a  structure  can  be  measured 
-  e.g.,  weight,  cost,  cost-effectiveness,  reliability.  While  lip-service 


has  been  paid  frequently  to  some  of  the  more  esoteric  merit  criteria,  the 
overwhelming  majority  of  the  work  performed  has  fallen  back  on  the  use  of 
weight  for  the  objective  function.  Weight  is  usually  important  in  most 
structures,  especially  aerospace  applications,  since  it  is  simple  to  calcu¬ 
late  and  it  is  noncontrover tible.  This  is  not  to  say  that  the  use  of  the 
other  criteria  has  not  been  investigated,  but  the  enormity  of  the  other 
aspects  of  the  structural  optimization  problem  has  generally  forced  re¬ 
searchers  to  choose  the  simplest  merit  criterion  initially  (and  finally) 
while  justifying  this  situation,  that  other  criteria  may  be  usable  within 
a  general  formulation  -  if  they  can  be  suitably  quantified.  The  use  of 
reliability  as  a  criterion  was  addressed  by  Moses,  Ref.  3,  where  concerns 
were  being  expressed  that  an  optimized  structure  is  more  likely  to  fail 
catastrophically  than  a  conventionally  designed  one.  Regrettably,  there 
appears  to  have  been  little  follow-up  work  along  the  lines  of  that  reference. 

If,  in  spite  of  the  expenditure  of  a  large  level  of  continuing  effort, 
major  fundamental  problems  still  exist,  what  is  the  nature  of  such  problems? 
Paradoxically,  it  appears  that  element  stresses  present  greater  problems 
than  deflections. 

Using  a  Lagrangian  formulation,  an  exact  solution  may  be  determined 
for  an  optimum  redundant  structure  subject  to  a  single  displacement  con¬ 
straint.  For  multiple  displacement  constraints,  the  same  approach  may  be 
used,  but  the  resulting  equations  become  nonlinear.  The  solution  of 
these  equations  has  proved  to  be  extremely  difficult  for  anything  but 
small-scale  test  problems.  For  stresses,  attempts  to  use  the  Lagrangian 
formulation  were  not  so  successful.  If  a  similar  approach  to  that  used 
for  deflections  is  applied,  the  Lagrangian  leads  to  fully  stressed  design. 
This  is  correct  for  determinate  structures  but  is  questionable  for  redundant 
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systems.  To  overcome  this  problem  the  effect  of  internal  force  distri¬ 
bution  with  variation  in  elastic  properties  must  be  introduced.  This  has 
been  accomplished  by  using  the  force  method  of  structural  analysis,  whereby 
the  redundant  structures  are  transformed  to  pseudo-determinate  systems  and 
the  internal  self-equilibrating  force  systems  are  treated  as  external  loads. 

The  basic  philosophy  of  the  force  method  of  structural  analysis  is 
well  rooted  in  antiquity,  definitely'  predating  the  much  more  popular 
displacement  method.  The  concepts  of  "cutting”  redundant  structures, 
applying  self-equilibrating  forces  and  computing  influence  coefficients 
have  been  universally  used  to  introduce  novice  structural  engineers  to 
the  mysteries  of  analyzing  complex  structures.  The  fundamentals  of  the 
force  method  are  addressed  in  many  publications  (see,  for  example.  Refs. 

4  and  5) . 

Although  force  method  concepts  are  particularly  simple  and  amenable 
to  a  ready  physical  comprehension  by  the  neophyte,  it  is  recognized  that 
some  computational  complexity  may  occur  in  poorly  formulated  problems. 

While  the  displacement  method  may  be  more  difficult  to  comprehend  initially 
(the  concept  of  a  stiffness,  per  se,  is  not  immediately  obvious),  the 
ability  to  automate  the  procedure  coupled  with  a  relative  lack  of  condi¬ 
tioning  problems,  has  led  to  its  dominance  in  the  field  of  structural 
analysis. 

In  the  early  work  in  structural  optimization  the  use  of  the  displace¬ 
ment  method  for  the  structural  analysis  stage  was  used,  more  or  less  auto¬ 
matically,  because  of  the  ready  availability  of  developed  finite  element 
programs  with  extensive  element  libraries.  Effort  was  concentrated  on 
the  development  of  rational  redesign  methods  using  various  strategies. 


Unfortunately,  the  most  successful  approach  via  the  Lagrangian 
formulation,  which  handles  displacement  constraints  in  a  very  efficient 

manner,  tends  to  break  down  when  handling  element  stress  constraints.  The 

S. 

problem  lies  in  the  expression  of  the  stress  constraint  in  the  form  — — — 1 

Ai  i 

where  is  an  element  force,  A^  its  area  or  thickness,  and  0.  the  al¬ 
lowable  stress.  In  the  Lagrangian  formulation,  derivatives  of  the  con¬ 
straints  with  respect  to  the  primary  variables  (A^)  must  be  formed.  Since 
in  a  redundant  structure  the  element  force  is  a  function  of  every 
element  area  A^ ,  a  term  8S^/8Aj  does  exist.  This  term  may  be  small,  but 
its  neglect  is  the  reason  the  Lagrangian  formulation  leads  directly  to 
fully  stressed  design.  In  the  displacement  method  no  analytic  expression 
for  8S^/9Aj  can  be  generated.  The  best  that  is  possible  lies  in  the  use 
of  finite  difference  approximations.  Since  the  optimality  criteria  form 
a  set  of  nonlinear  equations  which  can  only  be  solved  using  gradient  methods 
(e.g.,  Newton-Raphson) ,  the  generation  of  second  differentials  using  finite 
difference  methods  becomes  computationally  impractical.  The  solution  lies 
in  the  use  of  the  force  method  formulation.  In  this  approach,  the  redundant 
structure  is  represented  as  a  pseudo-determinate  one  in  which  S^  is  not  a 
direct  function  of  all  A  ^  .  The  true  variation  S^  is  reflected  by  including 
the  redundant  forces  as  variables.  By  this  artifice,  the  missing  term 
8S^/9Aj  is  incorporated  through  a  variable  redundant  force  term  which  shows 
up  as  a  compatibility  (displacement  restraint)  term  in  the  basic  Lagrangian. 
The  additional  terms  are  explicitly  differentiable  and  permit  the  generation 
of  nonlinear  equations  in  an  analytic  form.  These  equations  are  then  amenable 
to  solution  using  a  combination  of  linear-programming  and  Newton-Raphson 
techniques . 

Based  upon  this  force  method  approach  a  pilot  program,  OPTFORCE  I 
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was  developed  at  Bell  Aerospace  Textron  and  is  presented  in  Reference  4. 
Typical  of  pilot  programs,  OPTFORCE  I  was  limited  in  a  number  of  important 
respects.  It  was  one  of  the  principal  objectives  of  the  research  reported 
herein  to  expand  OPTFORCE  I  into  a  usable  computer  code  by  eliminating  or 
reducing  its  deficiencies;  the  new  code  is  aptly  named  OPTFORCE  II. 

Two  principal  options  were  available  to  accomplish  this:  Expand  the 
existing  OPTFORCE  I  analysis  program-  through  the  development  of  new  elements 
and  other  sophistications  or  acquire  an  existing  developed  program  and 
integrate  it  into  the  optimization  stage  of  OPTFORCE  I.  One  of  the  concerns 
in  acquiring  developed  programs  is  the  capability  and  efficiency  of  the 
available  programs.  Their  capacities,  as  measured  in  terms  of  element 
numbers  or  degrees  of  freedom  are  generally  adequate.  The  major  deficiency 
discovered  resided  in  the  area  of  frequency  constraints  expressed  in  the 
force  method  of  formulation.  Thus,  frequency  constraint  formulations  were 
developed  and  are  reported  in  this  text  as  well  as  additional  finite  element 
formulations. 

The  typical  form  of  a  frequency  constraint  is  that  the  basic  structural 
frequency  (undamped)  shall  be  greater  than  a  given  value.  In  some  instances 
a  less  than  constraint  may  be  specified,  although  this  is  less  likely. 

Higher  order  harmonics  are  generally  left  unconstrained.  While  a  frequency 
is  a  generalized  structural  response,  like  a  deflection,  its  incorporation 
into  an  optimization  program  will  introduce  additional  problems.  Their 
exclusion  from  standard  optimization  programs  may  be  indicative  of  both 
their  difficult  nature  and  of  their  lack  of  wide-spread  usefulness. 

Optimization  involving  frequency  constraints  of  various  types  has  been 
explored  (Refs.  5  &  6)  with  some  success.  Generally,  they  have  been  used 
alone  rather  than  in  combination  with  the  more  usual  stress,  deflection  and 


fabrication  constraints.  The  expanded  version  of  OPTFORCE  I  provides  an 
integration  of  frequency  constraints  with  the  other  types  mentioned.  A 
full  integration  of  the  frequency  constraints  into  an  optimization  program 
is  theoretically  similar  to  the  rigorous  treatment  of  any  other  types  of 
multiple  constraints,  i.e.,  the  problem  becomes  extremely  nonlinear,  but 
solvable  through  the  use  of  some  form  of  Newton-Raphson  technique. 

The  optimized  structure  may  be  more  sensitive  to  damage.  Hence,  it  is 
essential  to  make  provisions  for  the  service  reliability  or  vulnerability  of 
optimized  structures.  To  assess  the  service  reliability  of  the  structures, 
the  concept  of  partial  degradation  was  introduced.  This  may  simulate  battle 
damage  or  any  other  type  of  partial  failure  which  may  occur  as  a  result  of 
extensive  service  operation.  To  assess  the  suitability/reliability  of  the 
initial  structure,  the  behavior  of  the  degenerate  (damaged)  structure  is 
analyzed.  If  the  primary  structure  is  still  largely  capable  of  carrying  out 
its  prescribed  mission  even  when  damaged  then  its  reliability  is  acceptable. 
The  above  evaluations  of  reliability  are  clearly  only  qualitative  in  nature. 
Another  objective  of  the  study  program  was  to  establish  a  quantitative  means 
for  the  assessment  of  a  damaged  structure. 

To  determine  the  effect  of  damage  on  the  optimized  structure,  repeated 
re-analysis  of  degenerate  forms  of  the  original  highly  efficient  structure 
are  required.  Re-analysis  methods  have  been  under  development  since  the 
work  reported  in  1969  by  R.  J.  Melosh,  Ref.  7,  and  more  recently  with 
Increased  Intensity  by  other  researchers  such  as  V.  Venkayya  and  N.  Khot , 
Refs.  869.  The  search  for  more  efficient  re-analysis  methods  is  para¬ 
mount  to  the  solution  of  aerospace  structural  optimization  and  damage 
assessment  problems.  These  methods  are  vital  to  the  determination  of  the 
vulnerability  of  a  damaged  structure.  At  the  present  time  re-analysis 
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methods  can  be  broadly  classified  into  two  groups;  namely,  direct  and 
iterative.  The  former  eliminates  problems  associated  with  convergence  but 
could  be  burdened  more  with  problems  of  efficiency  than  the  latter.  Con¬ 
siderations  of  efficiency  are  of  utmost  importance  if  new  design  tools 
such  as  produced  herein  are  to  become  practical.  Examination  of  the  re¬ 
analysis  methods  developed  to  date  showed  that  viable  methods  do  exist  and 
some  are  applicable  to  large  scale  structural  analyses  and  design  of  air¬ 
frame  components.  These,  however,  were  not  cast  within  the  framework  of 
the  force  method.  Thus  a  new  method  of  rapid  re-analysis  of  damaged 
structures  was  developed  including  vulnerability  assessment.  Vulnerability 
methodology  must  account  for  the  interaction  of  optiraizat ion-reanalysis 
programs  in  order  to  provide  sufficient  details  of  the  damaged  structure 
such  that,  for  example,  residual  strength  can  be  determined.  Another  item 
of  vulnerability  is  the  change  in  the  vibration  characteristics  of  the 
damaged  structure.  Perhaps  what  is  of  most  importance  is  the  ability  to 
define  the  level  of  vulnerability  or  in  other  words  the  ability  to  deter¬ 
mine  the  operational  capability  of  a  damaged  structure.  This  facet  of 
vulnerability  was  Investigated  and  is  provided. 

As  a  result  of  the  efforts  previously  described  a  set  of  analytical 
tools  capable  of  designing  lighter  and  less  costly  structural  components 
is  available  to  the  designer.  It  was  appropriate  then  that  efforts  were 
expended  to  analyze  representative  structures  in  order  to  demonstrate  the 
utility  of  the  computer  codes  developed.  This  was  accomplished  on 
structures  varying  from  relatively  simple  trusses  to  a  swept  wlngbox  air¬ 
craft  structure.  Details  of  the  analyses  conducted  are  presented. 

The  analytical  tools  developed  and  based  on  the  force  method  of 
finite  analysis  increases  the  versatiblity  and  scope  of  analyses  which  can 
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be  performed  using  Che  force  method  of  approach  provided  in  OPTFORCE  I , 
Ref.  4.  However,  two  fundamental  questions  remain:  Does  the  force  method 
offer  an  improvement  in  efficiency  resulting  in  lower  computer  costs 
compared  with  the  displacement  method  of  finite  element  analysis?  And 
which  of  the  two  methods  is  the  most  accurate,  i.e.,  which  one  provides 
the  "correct"  minimum  weight  solution? 

As  noted  in  the  optimization  literature  when  the  iteration  algorithm 
is  based  on  the  displacement  method  of  finite  element  analysis,  a  complete 
re-analysis  of  the  structure  must  be  conducted  whenever  large  changes  are 
made  in  the  design  variables.  Thus,  substantial  computational  effort  is 
devoted  to  the  analysis  of  the  structure.  In  the  force  method,  the  basic 
framework  of  computational  effort  is  not  increased  with  changes  in  the 
design  variables  after  the  basic  and  redundant  load  system  is  selected. 
This  should  mean  that  there  is  an  improvement  in  efficiency  and  a  speedup 
in  the  iterative  analysis  for  the  force  method  algorithm.  Hence,  a 
comparison  of  computer  cpu  time  expended  by  the  OPTFORCE  II  program  solv¬ 
ing  classical  problems  with  known  solutions,  versus  an  existing  displace¬ 
ment  method  optimization  program  OPTIM  III,  Ref.  10,  solving  the  same 
problems  was  conducted.  This  provided  the  needed  comparison  data  to 
answer  the  fundamental  "dollars  and  cents"  question.  More  importantly, 
these  comparison  studies  also  provided  a  means  to  assess  the  relative 
accuracy  of  the  two  optimization  codes.  These  items  formed  the  basis 
for  recommending  future  development  of  a  general  purpose  optimization 
program. 

Thorough  discussions  of  the  technical  areas  discussed  above  are 
given  in  Sections  2.0  and  3.0.  Section  2.0  presents  the  theoretical 
background  of  the  OPTFORCE  II  code  and  emphasizes  the  expansion  of  the 
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pilot  program  OPTFORCE  1.  The  discussion  commences  with  a  review  of  the 
fundamentals  of  the  force  method  in  Section  2.1.  Derivation  of  finite 
element  matrices  for  bar  (axial  force),  membrane  triangle,  membrane 
quadrilateral  and  shear  panel  elements  within  the  context  of  the  force 
method  are  presented  next.  Weight  optimization  methodology  is  given  in 
Section  2.3  including  provisions  for  handling  variable  stress,  multiple 
displacement,  maximum  and  minimum  Size,  multiple  loading  and  natural 
frequency  constraints.  Rapid  re-analysis  and  damage  assessment  technology 
is  presented  and  illustrated  in  Section  2.4.  The  newly  developed  OPTFORCE  II 
code  Is  amply  exercised  in  Section  3.0.  Efficiency  studies  and  applica¬ 
tions  of  the  code  to  a  swept  wingbox  structure  are  described.  Program 
results  are  profusely  given  to  enable  the  reader  to  fully  understand  its 
capabilities.  It  is  noted  here  that  Volume  2  (Ref.  11)  of  this  report 
presents  details  of  the  program  and  the  input/output  format  of  OPTFORCE  II. 

A  technical  discussion  of  the  research  conducted  is  presented  in  Section 
4.0  wherein  conclusions  and  recommendations  are  given. 
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2.0  THEORETICAL  DEVELOPMENT 


1 


2.1  Fundamentals  of  the  Force  Method 

The  force  method  of  analysis  was  used  to  model  structural  behavior 
and  determine  minimum  weight  under  static  loading.  An  additional  requirement 
above  the  static  constraints  of  minimum  size,  stress  and  displacement  was  to 
constrain  the  structure  from  having  certain  of  its  fundamental  vibration 
frequencies  falling  within  a  specified  range.  In  order  to  accomplish  the  task 
of  stating  these  frequency  constraints  in  a  particular  mathematical  form 
compatible  with  the  weight  optimization  method  as  set  forth  in  Section  2.3, 
below,  it  was  necessary  to  formulate  a  force  method  dynamic  analysis  that 
would  yield  frequencies  and  inodes  when  the  structure  was  in  a  state  of  free 
vibration.  What  follows  is  a  discussion  of  the  theoretical  basis  for  both 
static  and  dynamic  analyses. 

2.1.1  Basic  Static  Analysis 

A  force  method  formulation  is  based  on  minimizing  the  complimentary 
energy  functional  with  respect  to  parameters  used  to  describe  a  chosen  approxi¬ 
mation  to  the  stress  field.  The  requirements  for  the  approximate  stress 
field  are  that  it  satisfies  equilibrium  throughout  the  body,  that  the  stress 
boundary  conditions  are  satisfied,  and  that  the  choice  admits  non-trivial 
solutions,  at  least  a  priori.  The  principle  of  minimum  complimentary  energy 
may  be  written  as 

6U*  +  6V*  -  0  •  (1) 


II 


U*  Is  Che  coaplimencary  strain  energy  and  V*,  for  a  finite 
element  idealization  of  a  body,  is  usually 

V*  -  -  {P}T  U)  (2) 

where  (P)  and  {A}  are  the  finite  element  nodal  forces  and  dis¬ 
placements,  respectively. 

The  stress  field  for  each  element  in  the  structure  is  given  by 
the  following  expression, 

«  [Nt]  {St)  (3) 

where  the  [n^]  are  functions  of  the  spatial  coordinates  defined  for  the 
element  and  {s^}  are  the  undetermined  parameters  for  the  stress  field, 
which  serve  as  coefficients  for  linear  combinations  of  the  [N^].  The  [N^] 
must  satisfy  che  homogeneous  equations  of  equilibrium,  and  must  admit  a  non¬ 
trivial  solution  for  each  stress  field  in  The  number  of  parameters  {S^} 

must  be  at  least  the  number  of  degrees  of  freedom  for  the  element  minus 
the  minimum  number  of  support  forces  needed  to  suppress  rigid  body  motion 
of  the  element.  With  the  constitutive  law  written  as 

Ej  *  [Ej]*1  0±  (4) 

the  complimentary  strain  energy  for  the  ith  element  is 

Ui  "  I  {Si}T  tfi3  tSi}  (5) 
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where  the  element  flexibility  matrix  is 

ffi]  -  fEj"1  fNtl  dV.  (6) 

The  total  complementary  strain  energy  of  the  structure  is  the  sum  of  the 
individual  strain  energies  of  the  elements.  If  (s)  is  the  collection  of  all 
the  {s^ }  and  [f]  is  a  matrix  formed  by  placing  the  [f^  along  the  main 
"diagonal"  of  an  otherwise  trivial  [f],  then  this  sum  may  be  written  as 

u*  -  \  {s}T  [fl  {S}.  (7) 

A  relation  between  {s}  and  {p}  must  be  found  in  order  to  use  (1) 
find  a  solution.  The  derivation  of  such  a  relation  presents  the  greatest 
difficulty  in  the  force  method  formulation  particularly  when  the  {s}  are 
devoid  of  physical  meaning.  One  approach  is  to  relate  the  {S^}  of  an  element 
with  {F^},  the  element  nodal  forces.  The  criterion  selected  for  this  relation 
is  that  {F^}  be  chosen  such  that  the  work  done  by  the  traction  field,  as 
calculated  from  the  stress  field  evaluated  on  the  element's  boundary,  is  the 
same  as  the  work  done  by  {F^}  through  the  nodal  displacements.  To  accomplish 
this,  it  becomes  necessary  to  choose  a  boundary  displacement  field  u^  in 
terms  of  the  nodal  displacements  {Aj};  thus, 

\  (8) 

The  work-equivalency  criterion  is  expressed  as 

{A1}T(rfY11T  [Lj  dS)  {St}  -{Aj}1  {  }  (9) 

s 

where  [  Lj]  is  [Nj]  evaluated  on  S,  Equation  (9)  yields  the  relation 

{^t}  *  [BJ  {St>  (10) 
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where 


r^i  -  ^  rY±7T  [l±]  as.  (id 

s 

The  Introduction  of  [Y^]  into  the  formulation  gives  it  a  hybrid  character 
(Ref.  12).  While  the  choice  of  u^,  appears  arbitrary,  an  interelement  continuous 
shape  would  be  best  to  choose. 

The  collection  of  equations  (10)  may  be  written  as 

{f}-[b]{S>.  (12) 

It  is  also  convenient  to  include  in  both  {F}  and  {s}  the  set  of  reactions 
{r}.  Equilibrium  equations  are  now  written  for  each  degree  of  freedom  as 

{p}-[cim  (i3) 

where  [c]  is  a  Boolean  matrix.  Substituting  (12)  into  (13)  yields 

{P}  -  (A]  [S).  (14) 

If  the  structure  is  statically  determinate,  then  [  A]  is  square  and  invert¬ 
ible,  and  (14)  is  solved  directly  for  {S},  and  thence  the  stress  field. 

In  order  to  find  the  nodal  displacements  {&},  equation  (14)  is  sub¬ 
stituted  in  (1)  and  variations  are  taken  with  respect  to  (P),  yielding  the 
relation 

(([AfV  [f]  [A]’1)  (PJ  -  {a}  .  (15) 

The  matrix  product  pre-mult ip lying  {P}  in  (15)  is  the  global  flexibility 
matrix  of  the  structure. 

If  the  structure  is  statically  indeterminate,  the  number  of  columns 
of  [A]  will  exceed  the  number  of  rows;  thus,  partition  (14)  into 

<w 


14 


where  [Aq]  is  square  and  Invertible.  The  subset  (x)  of  { S }  are  known  as 
the  "redundants".  Solving  (16  )  for  {Sq}  yields 

<so}  -  -[ao]_1  Ca13  {X}  +  UJT1  .  (17) 

Combining  (17 )  with  identity  relations  for  (x)  gives 

(S}  -  [bL]  (X>  +  [D]  {P}  (18) 

where 

[bj]  -  r-[Aorl  ;  [D]  -  '  (19) 

Equations  (18)  are  used  in  (17)  ,  and  then  substituted  into  (1)  ,  with 
variations  taken  with  respect  to  both  {x}  and  {p}.  The  two  sets  of  equations 
formed  are 

M  {X}  +  [*]  {P>  -  {0}  (20a) 

0]T  {x}  +  [a]  (P)  -  {A}  (20b) 

W  *  [b,]T  [f]  [bx];  M  *  [bj  [f]  [d];  [n]  =  [d]t  [f]  [d].  (2i) 

Equations  (20a)  are  solved  for  (x),  and  may  be  used  in  (18)  to  sub¬ 
sequently  determine  the  stress  field  throughout  the  structure.  Substitut¬ 
ing  (20a)  for  {x}  in  (20b)  yields  the  global  flexibility  relations 

OH  {P}  -  (A)  (22) 

where  the  global  flexibility  matrix  is 

OH  -  M  -  l*f  M”1  M.  (23) 

Note  that  since  {P}  is  known  and  (A)  is  unknown  that  (15)  or  (22)  need 


not  be  inverted;  on  the  other  hand,  the  inversion  of  (Aj  requires  the  same 
amount  of  effort  as  the  inversion  of  a  global  stiffness  matrix. 

In  order  to  facilitate  incorporation  of  the  above  formulation  into  the 
weight  optimization  method,  several  special  properties  of  the  force  method 
need  to  be  recognized.  First,  the  finite  elements  used  have  volumes  which 
can  be  written  in  the  form 


Vi  “  Ai  Ri 


(24) 


where  R^  is  a  one  or  two-dimensional  element  domain  and  A^  is  the  correspond¬ 
ing  two  or  one-dimensional  element  size.  For  a  given  structural  layout  and 
discretization,  will  be  fixed  for  an  element  with  A^  to  be  determined. 
Furthermore,  the  bounding  surface  of  the  element  can  be  written  as 


-  At  Ct  (K) 

where  is  a  one  or  zero  dimensional  ’’curve"  surrounding  R^. 

Suppose  that  the  parameters  for  the  stress  field  {S^}  are  dimensionally 
given  in  units  of  force.  This  seems  logical  if  the  set  of  reactions  are  to 
be  incorporated  with  {S};  furthermore,  matrices  fB^l  and  thus  [A]  become 
strictly  non-dimensional  under  this  assumption.  If  shape  functions  (N^J 
are  defined  as 

t  N±1  -  At  [Nj  (26) 

then  use  of  (26)  and  (24)  in  (6)  yields 
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<fil 


/  (t“  t N.J T)  (El”1 

Rj  i 


*71S1»  *: 


dR. 


(27) 


where  IfJ  is  Independent  of  A^.  In  addition,  (26)  and  (25)  are  used 
in  (11)  to  form 
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(28) 


*B±1 


f  ,YilT<»lfLil>  Ai  dci 


where  (L^l  is  l N^]  evaluated  on  C^,  thus  rendering  (B^  independent  of  A.. 
The  result  is  that  (A],  and  thus  [b^]  and  [  D]  are  independent  of  the  element 
sizes  and  strictly  a  function  of  the  geometric  layout  and  discretization 
of  the  structure.  This  eliminates  repetitive  computations  as  element  sizes 
are  changed  during  the  weight  optimization  analysis. 

The  second  special  property  of  the  formulation  becomes  apparent  when 
equations  (18)  are  partitioned  by  elements;  thus, 

{S1)  =  [b1  ]  {X}  +  {Dj  {P}  .  (29) 


In  addition,  the  reactions  are  calculated  as 


{R}  -  [b1  ]  {X}  +  I Dr1  (P)  .  (30) 

R 

The  diagonal  nature  of  ff],  combined  with  the  partition  (29)  and  (27) 
yields  for  the  expressions  (21) 


[♦1- 


£ 

i* 


1  i 


£  - 
i=l 


Inl  -  z  [i!,] 

i-1  i 


(31) 


where  N  is  the  number  of  elements  and 
£ 

[♦t]  =  Cb1]T[fi][bi];[j;i]  -  -  [D1]T[f1][D1]  .  (32 

Note  that  [$  ],  [^]  and  are  independent  of  A^,  and  that  the  A^ 

dependence  of  the  basic  static  equations  (20)  may  be  explicitly  written. 
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2.1.2  Basic  Dynamic  Analysis 

As  stated  above,  a  force  method  approach  to  solutions  of  structural 
problems  requires  that  equilibrium  be  satisfied  throughout  the  body.  In 
order  to  extend  the  force  method  to  problems  of  structural  dynamics,  it 
is  required  that  the  appropriate  forma  of  Newton’s  second  law  be  satisfied 
throughout  the  body.  It  will  be  convenient  to  use  D'Alembert's  principle 
and  introduce  "inertial  loads"  into  the  equations  of  "dynamic  equilibrium". 

When  analyzing  a  structural  dynamics  problem  using  a  finite  element 
lumped  mass  approach,  the  structure  is  assumed  to  be  composed  of  particles 
at  the  node  points  containing  mass,  connected  by  massless  elements.  Since 
the  elements  will  have  no  inertial  loading  as  a  result  of  this  idealization, 
the  equations  of  dynamic  equilibrium  reduce  to  those  of  static  equilibrium. 
Thus,  the  static  shape  functions  [N£]  may  still  be  used  to  describe  the 
stress  field  within  the  element.  The  formulation  proceeds  exactly  as  in 
the  statics  case  from  equations  (3)  -  (12). 

Equations  (13)  must  be  modified  to  include  the  inertial  loading 
(Pj)  as  part  of  the  overall  "external  load".  With  {P}  designating  the 
actual  applied  loading  only,  (13)  is  changed  to 

{P>  +  (Pj)  -  [c]  {F>  .  (33) 

The  formulation  proceeds  much  like  the  statics  case,  except  {p}  +  {P£}  is 
substituted  for  {P}  in  the  equations  where  it  appears.  In  particular 
(22)  (or  15)  becomes 

[3]  ({P}  +  (Pj))  -  (A).  (34) 

The  Inertial  loads  are  related  to  the  displacements  by 

(Pj)  -  -[M£]4‘}  (35) 
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where  [m^]  is  the  diagonal  lumped  mass  matrix  and  a  dot  indicates  a  time 
derivative.  Integrating  (2.35)  with  respect  to  time  from  0  t n  t  gives 


(A(t)  )  -  {AoI-Cm^'V  {Pj (t ) }  dr 


(36) 


where  {A0 )  are  the  initial  nodal  velocities.  Integrating  again  yields 


(A(t) }  =  {A0}  +  {A0)t  -  [M£r10/t  (t-T)  {P  (T) }  dT*  (37) 

where  CA0 }  are  the  initial  nodal  displacements.  Substituting  (36)  ■  into 
(34)  and  rearranging  terms  yields 


13}  <Pj}  +  CMJir10/t(t-T)  {Pi(t)>  dr  =  {A0>  +  (A0 )  -  {Astat(t)l 

where 


{Astat(t)}  *  M  {P(t))  *  (39) 

The  basic  equations  (38)  are  integral  relations  with  (Pj)  as  the  basic 
variable.  Differential  equations  with  (A)  as  the  basic  variable  could  be 
derived,  but  that  type  of  formulation  would  not  be  in  the  spirit  of  the 
force  method,  where  forces  are  the  basic  unknowns.  Furthermore,  the  form 
of  (38)  is  directly  analagous  to  a  stiffness  formulation. 

Fundamental  modes  and  frequencies  for  the  structure  may  be  found 
by  formulating  the  free  vibration  problem.  Assume  that  (P^.}  may  be 
written  as 

{PT}  *  {P.  }  cos  tut  +  (PT  )  sin  u)t  .  (40) 

I  Ic  Is 

*o/tC0/t"{P1(T))dT]dt'-  t0'/C  {Px(T)}dT  I1-,/1  t'{PI(t"))dt' 

0 

-  t  ,c{p  (x)}dT  -0/Ct(Pt(t) }dT 

o'  1  1 

-  0/t(t-T)  {Pi(t)>  dT 


(38) 
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Using  (40)  in  (38)  ,  and  remembering  chat  for  free  vibration  (P(t)}  -  {o}, 
gives 

(Cfl  -  -7  {pi}  +  (T  K3"1  {PIc}  '  {A«}) 

U)  0) 

(41) 

+  (“  [M^]"1  {PIg}  -  (A0»t  -  (01  . 

With  (A0 }  and  (A0)  arbitrary,  {P^}  and  (P^)  may  be  chosen  such  that 
(41)  reduces  to 


(CM  -  — 2[M£]"1)  fPj>  =  {0}  .  (42) 

In  order  for  non-trivial  (Pj)  to  exist,  the  determinant  of  m  -  -  cm.]-1 

ur 

must  be  zero.  The  resulting  equation  yields  the  fundamental  frequencies. 

In  order  to  facilitate  the  numerical  calculation  of  modes  and  frequencies 
by  using  standard  matrix  operation  subroutine  packages,  it  is  desirable 
to  make  the  substitution 


{Pj}  -  [M£]Ss  {q}  (43) 

in  (42).  Note  that  [MgT*  i»  easy  to  calculate,  considering  the  diagonal 

L 

nature  of  [m^].  Pre-multiplying  (42)  by  [M^J5  gives 


where 


(Cq3  -  ~W\)  <q>  -  {0} 

0) 

[Q]  -  [m^CM  . 


(44) 


(45) 


Thus,  the  fundamental  frequencies  are  the  inverse  square  roots  of  the 
eigenvalues  of  [qJ. 
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The  eigenvalues  of  [Q]  are  found  simultaneously.  Suppose  the  kth 
eigenvector,  i.e.,  the  one  corresponding  to  frequency  w^,  is  denoted  as 
{qk> ;  then,  the  kth  inertial  load  mode  {P^}  is  found  from  (43)  ;  the  kth 
redundant  "force"  mode  from  (20a)  the  kth  parametric  mode  from  (18)  ;  and 

the  k  stress  mode  from  (3)  •  The  ktn  displacement  mode  is  found  most 
easily  by  comparing  (42)  with  (34)  to  yield 

{Ak>  -  [M^]"1  {P^}  .  (46) 

0) 

Mode  normalization  may  take  place  on  either  (P^  >  ~  fV  • 

The  diagonal  mass  matrix  [M^]  may  be  found  as  a  function  of  the 


element  sizes.  Introducing  the  mass  per  unit  volume,  p^,  the  mass  of  the 
i*"*1  element  is  written  as 


mi  *  pi  Vi  *  Pi  Ai  Ri  “  Ai  mi 


(47) 


where  use  has  been  made  of  (24)  ;  thus,  ti^is  known  from  input  data.  This 

mass  is  equally  divided  and  assumed  lumped  at  the  node  points  of  the 

element.  The  number  8  ,  is  the  fraction  of  the  mass  of  element  i  located 

gi 

at  node  g.  Note  that  Bg^  can  be  defined  as  zero  for  those  nodes  that  are 
not  part  of  element  i.  Arbitrary  criteria  for  assignment  of  to  the 
various  nodes  may  be  made.  In  the  present  study,  8  .  *  1/N  •  where  N 
is  the  number  of  nodes  associated  with  element  1,  was  chosen.  If  g'  is 
a  degree  of  freedom  associated  with  node  g,  then 


ne 

V  "  £  Bgi  Al 


(48) 


Problems  of  forced  vibration,  as  well  as  a  consistent  mass  formula¬ 
tion,  may  be  derived,  but  were  not  utilized  in  the  present  study.  The 
Interested  reader  is  referred  to  Reference  13. 
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Derivation  of  Element  Matrices 

Material  and  geometric  information  combined  with  an  admissible  choice 
for  a  stress  field  will  produce  matrices  for  a  given  element.  In  fact,  only 
two  matrices  are  of  importance  (at  the  element  level):  [f^]  and  [B  ].  Both 
of  these  are  defined  relative  to  local  coordinates.  The  [B^]  matrix  must  be 
transformed  to  global  coordinates  before  Incorporating  it  into  [b]  for  use 
in  equation  (2.12). 

2.2.1  Truss  (Rod)  Element 

Figure  1  shows  a  typical  truss  element.  The  x  axis  indicated  is 
a  local  coordinate  axis.  The  element's  domain  consists  of  the  points 
0<x<L,  while  the  "bounding  surface"  is  simply  the  node  points.  A  single 
parameter  S  is  used  to  describe  the  stress  field: 


Fi,4i 


Figure  1  Truss  (Rod)  Element 
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where  A  Is  Che  cross  sectional  area  of  Che  truss  element,  and  serves  as  the 
element  size.  Thus 

[N]  -  1  (50) 

wi  th 

[E]'1  »  |  •  (5D 

For  this  element. 
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$ 


(56) 


Figure  2  Plane  Stress  Triangle 


where  t  is  the  cross-sectional  thickness  of  the  element,  and  serves  as  the 
element  size.  A  is  the  triangular  surface  area,  whose  square  root  serves 
as  a  linear scaling  factor.  The  [N]  matrix  is  simply 

[53  -  —  [i]  .  (57) 

Ja 

Orthotropic  material  behavior  will  be  assumed.  When  the  material  axis, 

x  ,  is  aligned  with  the  x  axis,  the  3x3  compliance  matrix  is  written 
m  c 

in  the  form 
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(59) 


where 


v  v 

A  a  SL  « 

P  E  E 

x  y 


Use  of  (59)  and  (57)  in  the  calculation  of  [f]  yields 


[f]  -  [E]  .  (61) 

The  calculation  of  [b.]  is  performed  by  summing  the  contribution 

1. 

from  each  edge  of  the  element.  Figure  3  represents  a  typical  edge  of  a 
triangular  (or  any  general  polygonal)  element  with  nodes  (i,j)  at  the  end 
of  the  edge.  The  displacement  fields  u^_^  are  defined  to  be  linear  in  the 
edge  coordinate  x' ,  thus 


Vj  <x'> 

Vi  <x'> 


0  0 - 1  -  J-  0  ~  0 - 0  0 


0  0- 


1  -  X 1  0  2Sl  no 

JT  0 


where 


T 

(A)  ■  [A  ,A  •••A  ,A  ...A  ,A  ***A  ,A  J, 
X1  ^1  Xi  ^i  xj  ^ 


Figure  3  Edge  Coordinate 
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N  is  the  number  of  nodes  for  the  polygon  element.  Note  that  [Y^]  is  a 
(2  x  2N)  matrix,  and  thus  for  the  triangular  element  it  would  be  of  order 
(2  x  6). 

To  determine  the  traction  field  for  an  edge,  equilibrium  must  be 
satisfied  for  an  edge  segment.  Noting  Figure  4  ,  the  traction  field  is 
determined  as 


Tx  (x’) 
Ty  (*') 


sin  9  0  -  cos  0 

0  -  cos  0  sin  0 


10  dx*  cos  9 

r  y 


T  dx'  sin© 
XU  I 


O  dx'  sin  0 
x 


T  dx’  cos  0 


T  dx’ 
x 


Figure  4  Edge  Segment  Free  Body  Diagram 


where  the  stresses  are  evaluated  on  edge  i-j .  For  the  constant  stress 
triangle,  the  relation  is  (2.56);  thus. 


[l*3  *  ~ 


1  sin 

Ja  0 


0  -cos  0 

-  cos  0  sin  0 


Using  (64)  and  (62)  in  one  term  of  (53)  ,  say,  edge  1-2,  gives 


-viaiKsc 
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-  f 


i-x*/i'  o 

0  l-x'/A' 

x’/A’  0 

o  x'/r 

0  0 

0  0 


1 

/X 


sin  6  0  -  cos  9 

0  -  cos  9  sin  9 


dx’ 


£  ®l-2^  2<Jk 


y2-yi 

0 

-  (x2~x1) 

0 

-  (X  -x  ) 

y2'yl 

y2-yl 

0 

-  (x2-x.) 

0 

-  (Xj-Xj) 

y2"yl 

0 

0 


0 

0 


0 

0 


(65a) 


where  the  subscripts  pertain  to  the  local  coordinates  of  the  elements  grid 
points. 


Similarly, 


0 

0 

y3-y2 

o 

y3-y2 

o 


0 

0 

0 

0 

0 

-  (x3~x2) 

(x3-x2) 

y3-y2 

0 

-  (x3-x2) 

/■> 

N 

X 

X 

w 

1 

fo 

yry3 

o 

o 

0 

yry3 


-  (Xj-Xj) 

0 

0 

0 

-  (Xj-Xj) 


“  (Xj-Xj) 

yry3 

0 

0 

-  (Xj-Xj) 

yl“y3 


(65b) 


(65c) 
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Adding  (65a-c)  to  form  C B^]  gives 


t**1  ■  i, 


2<Jk 


y2-y3 


y3~yl 


yry2 


V*2 


xl~x3 


X2"X1 


X3'X2 


y2-y3 


Xl"X3 


y3-yl 


X2~X1 


yry2 


o 


2.2.3  Plane  Stress  Quadrilateral 

Figure  5  shows  a  typical  quadrilateral  shaped  plane  stress 
element.  The  x  and  y  axes  shown  are  considered  to  be  the  elements  local 
axes.  A  minimum  of  five  parameters  are  necessary  to  describe  the  stress 


Figure  5  Plane  Stress  Quadrilateral 


field.  The  admissible  functions  chosen  were 


(67) 


T  • 

where  { S}  »  S2>  S^,  S^,  S^]  and  [K]  is  the  (3  x  5)  matrix  on  the 

right-hand  side  of  (67).  The  compliance  matrix  [E]  1  is  the  same  as 
that  given  for  the  plane  stress  triangle  element.  Performing  the  calcula¬ 
tion  for  [f]  yields 


[f] 


(E"1) 


11 


y"  (E'1) 


/A 


11 


(E-1) 
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/A 


C  (E-1) 


21 


(E"1) 


31 


xx 

.2 


(E_1) 21 


—  symmetric — 


(E'1) 


>Ta 


22 


«\l  /  (E_1)22  ^  (E"1)22 


h  CE-1) 3l  (E-1)32 
vA 


Xc  (E~1)v>  ._-l. 

jx  32  <E 


(68) 


where  A  is  the  surface  area  of  the  quadrilateral;  x^  and  are  the 

coordinates  of  the  centroidal  position;  and  I  ,  I  and  I  are  second 

xx  yy  xy 

moments  of  area. 
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As  with  Che  plane  stress  triangle,  the  edge  displacements  are 
assumed  linear.  In  addition,  the  relation  (63)  holds,  but  now  the 
evaluation  of  the  stresses  are  more  complex  due  to  the  coordinate 
dependence  as  specified  in  (67).  In  order  to  express,  for  example, 

the  stresses  on  edge  1-2  due  to  x*  dependence,  it  is  worth  noting  that 
the  relation  between  x  or  y  and  x'  is  linear  and  thus  can  be  specified 
much  like  the  displacement  along  the  edge: 


rn  x*  -i  x' 

x  “  X1  [1  “  IT3  +  *2  V 

y  =  Yi  [1  -  frl  +  y2  fr 

Thus  for  edge  1-2 


(69) 


BmI-  L 


sin0  0  -cos© 


-cos0  sin0 


-i—  ify  (i_  ~.)+y  0  0  0 

A^V1  *,;“2  fc,J 

0  0  J-  ^Cxj  (1-  p-)+x2jr]  0 


1 

/A 


~  sin0  ~  sin0  [y^l-  f^)+y2  frf  0 


/A 


/a 


:os0 


Vr030"  koaQ[xi(1~  i')+*2  r']^ine 


(70) 
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2.2.4  Symmetric  Shear  Panel 


The  thoretical  formulation  for  a  shear  panel  may  be  approached  in 
various  ways.  One  way  would  be  to  consider  it  to  be  a  special  case  of  a 
plane  stress  quadrilateral  with  normal  stresses  identically  zero.  Thus, 
for  a  general  quadrilateral  shear  panel 


T 


xy 


(74a) 


[f]  =  (E-1)33 


(74b) 


C*il- 


(74c) 


Equations  (74  a-c)  are  taken  from  (67),  (68),  and  (73)  where 

(here  identified  as  S)  is  the  only  non-zero  parameter.  The  element 
apparently  has  eight  degrees  of  freedom,  but  only  four  are  independent; 
or,  shifting  the  argument  to  a  force  method  point  of  view,  the  loads 
on  the  panel  must  be  such  that  the  constraints  implied  by  (74c)  are 
satisfied;  that  is,  a  panel  that  is  supported  against  rigid  body  motion, 
may  not  be  loaded  in  an  arbitrary  manner  at  the  free  nodes  and  still 
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maintain  its  Identity  as  a  shear  panel.  It  is  quite  possible  that  a 
structure  assumed  to  be  composed  of  quadrilateral  shear  panels  would 
yield  equations  (14)  with  dependent  or  contradictory  equations. 

Since  the  primary  application  of  quadrilateral  shear  panels  is 
modeling  spar  webs  in  wing  structures,  the  most  obvious  solution  to  the 
problem  described  above  is  the  creation  of  a  symmetric  web  element. 
Figure  6.  In  this  form,  only  half  an  element  is  used  in  the  analysis, 
and  thus  only  two  nodes  remain.  When  rigid  body  modes  are  restrained, 
only  one  degree  of  freedom  will  remain.  Note  that  some  support  condi¬ 
tion  on  nodes  1  or  2  implies  some  loading  on  the  lower,  non-analyzed 


Figure  6  Symmetric  Shear  Panel 


half  of  the  web.  The  element  matrices  may  now  be  written  for  the  (upper 
half  of  the)  sytanetrlc  shear  web  as 


T 

xy 


£  sXj-Xj) (yj+yj)]' 


(75a) 
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[f] 


1 

G 

xy 


Cb±3  **  H 


[(x2-x1/(y1+y2)]5s 

-[(y1+y2)^x2_xi^,S 

[(x2-x1)/(y1*y2)]is 

^(yl+y2)/<X2“Xl)^ 


where,  in  (75b)  the  x  axis  is  chosen  as  the  material  axis. 


(75b) 


(75c) 


2.3  Weight  Optimization  Method 

In  this  section,  the  optimization  method  developed  is  described.  What 
is  desired  is  to  obtain  the  minimum  weight  of  a  structure  of  given  geometric 
layout,  material  properties  and  loading  that  satisfies  a  set  of  size  and 
structural  constraints.  The  basic  unknowns  of  the  optimization  procedure 
are  the  element  sizes,  A^,  i»l,  ...,  N£,  where  is  the  number  of  elements. 
2.3.1  Theoretical  Basis 

The  theoretical  basis  for  the  generation  of  the  proper  equations 
for  the  optimization  procedure  are  the  Kuhn-Tucker  optimality  criteria 
(Reference  14).  These  criteria  are  formulated  by  extending  the  standard 
"Lagrange  multiplier"  method  of  optimizing  a  quantity  subject  to  constraints 
to  cases  of  inequality  constraints.  An  informal  review  of  the  theory  is 
now  presented. 

For  the  basic  Lagrange  scheme,  a  function  Z*f  (x.  ,  x,,  . ..,  x  ) 

i-  i  n 

subject  to  equality  constraints  g, (x  ,  x_,  ...,  x  )  ■  0,  i-1,  2,  ....  m<n, 

i  i  i  m 

is  to  be  optimized.  A  quantity  known  as  the  Lagrangian  is  formed  as 
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m 


1.  -  f  (x  )  +  EX  g.  (x  ) 

J  i-1  1  3 

Derivatives  are  taken  with  respect  to  each  x^  and  X^and  set  to  zero  to 
obtain  the  system  of  equations: 


(76) 


3L 

3x, 


11  +  r  X .  ffl 

3xj  1.-1  SXj 


j-1,2, 


(77a) 


3L 

ax. 


i-1, 2, . . . ,m 


(77b) 


Equations  (77b)  are  just  a  reiteration  of  the  constraint  equations.  The 
system  (77)  contains  (m+n)  equations  for  (m+n)  unknowns.  Note  that  there 
are  no  restrictions  on  the  signs  of  either  the  x's  or  the  X's. 

Suppose  the  optimization  problem  contained  inequality  constraints, 
say,  less  than  or  equal  to  constraints,  and  the  function  Z  is  to  be 
minimized  subject  to  those  constraints.  Then  a  set  of  variables  may 
be  introduced  to  make  the  inequality  constraints  equality  constraints  by 


g1+Gi  -0 


Now  the  Lagrangian  is 


i  —  1,2,... ,m. 


(78) 


l  ■  f  +  IX.  (g.  +  G.  ). 
i-1  111 


Taking  derivatives  with  respect  to  Xj's,  X^'s  an^  G^'s  gives 
3L  Jf  "  ^  n 

3xj“  i_1  Ai  3xj  "  U 


(79) 


(80a) 
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(80b) 


Si  +  Gi2  *  ° 


3Gi  "  2XiGi 


0  . 


(80c) 


Equation  (80c)  implies  that  either  X^  or  equals  zero.  If  G^  *  0, 
then  g^  -  0,  that  is,  the  constraint  is  satisfied  at  equality,  or,  as 
it  will  be  described  later,  the  constraint  is  "on".  If  X^»  0,  then 
would  not  necessarily  be  zero  and  the  corresponding  constraint  is  g^lO, 
or  satisfied  in  the  inequality  sense,  or  "off". 

Since  a  minimization  is  required  here,  the  pure  second 
derivatives  of  L  must  be  non-negative;  in  particular,  from  (80c). 


32L 


2^)0 


(81) 


or,  that  the  Lagrange  multipliers  must  be  non-negative.  The  discussion 
above  alL 
(80b,  c): 


above  allows  for  the  elimination  of  G^  from  the  analysis  by  writing  for 


If  X^  -  0,  then  g^  <0; 
If  X^  >  0,  then  gt  «  0; 
X^  non-negative. 


(82) 


Greater- than-or-equal- to  constraints  may  be  multiplied  by  -1  to  convert 
to  less-than  constraints,  or  their  corresponding  multipliers  must  be  non¬ 
positive.  Since  a  strict  equality  constraint  can  be  thought  of  as 
simultaneously  a  less-than  and  greater-than  constraint,  the  Lagrange 
multiplier  will  be  unrestricted  in  sign,  as  seen  above  in  (77). 
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It  is  necessary  for  the  generalized  problem,  then,  to  satisfy 
equations  (80a)  and  conditions  (82)  at  the  optimal  condition.  There 
may  be  other  extremum  where  these  conditions  are  satisfied  as  well.  Only 
when  certain  conditions  of  convexity  are  satisfied  are  these  conditions 
sufficient  for  an  optimum;  furthermore,  the  equations  themselves  give  no 
clue  on  how  to  solve  them. 

2.3.2  Problem  Formulation 

The  solution  method  for  a  system  of  equations  such  as  (80a)  and  (82) 
can  only  be  determined  after  the  specific  form  of  the  system  is  known.  The 
weight  optimization  problem  is  thus  formulated.  The  total  weight  of  the 
structure  must  be  minimized;  hence,  the  function  Z  is 

NE  _ 

M  =  ^  Wi  Ai  (83) 

where  is  the  weight  per  unit  element  size,  and  may  be  calculated  as 

Wj  *  n^g  (84) 

where  m^  is  defined  in  (47)  and  g  is  the  acceleration  due  to  gravity. 

In  addition  the  structure  is  subject  to  the  following  constraints: 

a)  Minimum  size  constraints:  The  element  sizes  (thickness,  cross- 
sectional  area)  must  exceed  a  certain  minimum  (positive)  value  due  to 

constraints  arising  from  manufacturing  capabilities.  If  the  minimum 

* 

value  for  the  ith  element  is  Aj,  then 

Ai  >  Ai  l**l»  2.  ....  Ng-  (85) 

b)  Maximum  Stress  Constraints:  Suppose  the  yield  (or  some  other 
failure)  stress  for  the  ith  element  is  o^.  Then  a  stress  constraint  must 
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be  in  the  form 


0  <  a*  i  «  1,  2,  ....  N  (86) 

ei  1  11 

where  o  is  an  effective  stress  measure.  Not  only  does  0  depend  on 
ei  ei 

the  element  sizes,  but  it  also  depends  on  the  state  of  loading. 


c)  Displacement  Constraints:  A  typical  displacement  constraint 


is  in  the  form 


4j  |  “j  J  ■  '■  2 . hdc  1 

*  th 

where  A^  is  the  limiting  value  on  the  j  displacement  that  is  desired 
to  be  constrained.  The  number  of  such  constrained  displacements  is  N^. 

d)  Frequency  Range  Constraints:  In  general,  these  constraints 
may  be  written  as 


^  t  -  1,  2,  ...»  ( 

til 

where  is  the  k^'th  fundamental  frequency  to  be  constrained  to  the  l 
*  th 

extreme  value  for  the  l  constraint  in  the  set  of  N^  frequency 
constraints. 

While  the  mathematical  forms  of  the  constraints  are  rather  simple 
to  write,  the  actual  calculation  of  the  left-hand  sides  in  terms  of  the 


element  sizes  are  not  simple  for  either  force  or  displacement  method 
formulations.  The  form  of  (80a)  demands  that  an  explicit  differentia¬ 
tion  of  the  constraint  equations  (when  stated  in  a  specified  form)  with 
respect  to  the  element  sizes  be  performed.  In  Section  2.1  of  this  volume, 
the  dependence  of  many  of  the  basic  matrices  on  the  element  sizes  is 
explicitly  given.  The  state  of  the  applied  loading,  of  course,  is  also 
given;  however,  the  set  of  redundants,  (x) ,  vhlch  acts  as  a  conduit  of 
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information  to  find  the  displacements  and  stresses,  cannot  be  written  in 
a  direct,  explicit  manner  as  a  function  of  the  element  sizes.  It  may  be 
advantageous  to  extend  the  set  of  basic  variables  to  include  { X } ;  and,  if 
done,  the  relation  between  (x)  and  element  sizes  must  be  explicitly  written. 
This  relation  is  nothing  more  than  (20a)  ,  and  will  serve  as  another  set 
of  constraints.  Note  that  the  explicit  element  size  dependence  in  (20a) 
is  one  which  is  a  sum  of  terms  each  containing  an  element  size  in  the 
inverse  first  power.  The  displacements  are  calculated  from  (20b)  ,  where 

this  condition  also  exists,  as  well  as  (3)  used  to  calculate  the  stress 
field.  It  thus  seems  likely  that  constraints  (86)  and  (87)  will  also 
be  in  this  form.  Equations  (85)  can  be  converted  into  this  form  as  well. 
Excluding  (88)  for  now,  the  constraint  equations  may  be  rewritten  in  one 
of  two  forms: 


nE  r  1 

e  ru  -  i  *  o 

1-1  Aj 


or 


NE  c  2 

£  _ii 


-  o  . 


j=l  A 


j 


Substituting  (89)  and  (83)  into  (76)  yields 


L 


NE  - 

E  V  A 
i-1  1  1 


ml  ! 
E  A1 

i-1 


NE  c  1  m2  2  Ne 

(  E  ^1-1)  +  E  X  (E 

j-1  Aj  i-1  j-1  A 


where  m^  and  are  the  number  of  constraints  of  the  type  of  (89a)  and 
(89b)  ,  respectively.  Differentiate  (90)  with  respect  to  A^  to  yield 


9L 

9Ak 


C 


2 

ik  ' 


(89a) 


(89b) 


(90) 


(91) 
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Now  form 


N 


E  3L  ?E  -  .  .  “l  ,1  ,  C  1 


m2  2  c  ^ 


L  +  E  \  H  “  1  W1  A.  +  E  X  <  E  -  1)  +  I  X*  E  _ii 
*  dAk  1=1  1  i-1  j-1  A  1=1  j=l  Aj 


N* 


*?E  l  ”l  ,2  C,1  -  ZE  £2  X?  Cik  -  2W  -  Z1  X? 


Sri  *  *•  «*» 


i-1 


1=1 


(92) 


Evaluating  (92)  at  the  optimum  condition  will  imply 


L’w*:  V  0,k"1,  *",Ne 


(93) 


where  W  is  the  optimum  weight;  thus,  substituting  (93)  into  (92)  gives 


I1  xj  -  w* 
i-1  1 


(94) 


at  the  optimum  design.  Thus,  the  sum  of  the  Lagrange  multipliers  for  the 
constraints  written  in  the  form  (89a)  equals  the  minimum  weight  of  the 
structure  for  the  optimum  design.  This  circumstance  will  play  a  prominent 
role  in  the  solution  algorithm. 

The  constraint  equations  are  now  rewritten  to  be  in  one  of  the 
forms  shown  by  (89) . 

a)  Minimum  Size  Constraints  -  Equations  (85)  are  rearranged  to 

form 

* 

g  «  \  ■  liu  L  =  1,  2,  ....  N  .  (95) 

A.  *  “ 
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M  stress  arna  -  *-  “ 

dependent  on  the  physical  — “**'*  f“  * 

a,  the  mathematical  for.  of  the  stress  field  in  terms  of  the^ehosen 
stress  parameters.  In  general,  .  Mlses-Bench,  criteria,  <»*  «WV 
Will  be  used  for  all  the  elements, 
i)  Trims  Element 


xy 


.  -<o^  -  -f  <^-xriM- 

■i  1  1 


(96) 


tl)  Plane  Stress  Trianfele 

2  .  .2 


(<Jx  +  V°x°y  +  3txy)1S 
ei 


2A 


-V  <»1s1-aiVjsr) 

■■iF  ■ 


(97) 


iii)  Plane  Stress  Quadrilateral 


2.  _  2 


+  3t  _.2) 


°e  *  <°x  +  °y  "  °x  °y  '  *y 


(98) 


Note  that  for  this  element, 


the  stress  is  evaluated  arbitrarily  at  the 


centroid  of  the  quadrilateral. 
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iv)  Shear  Panel 


(103) 


where 


“ij 


+  ^x  "i  NDN  ~i 

-  *»  '  +  £  "i*  P‘  • 


(104) 


d)  Frequency  Range  Constraints  -  In  order  to  formulate  these 
constraints  correctly,  it  is  necessary  to  rewrite  equations  (20)  for 
free  vibration.  The  k1'*'  fundamental  mode  and  frequency  satisfies 


[♦]  +  [*]  (Pk}  -  {0} 


WT  (V  +  <[03  -  A  Em*]"1)  {Pk>  -  {0}  . 


(105) 


T 

Premultiplying  the  first  equation  by  {j^}  ,  the  second  equation  by 


{Pfc}  ,  and  adding,  yields 


\}r  M  +  2  [*]  (pkJ  +  (Pk}T  [03  {Pk> 


^  {pk}T  ^  <Pk>  -  o  . 


uk 


(106) 


The  first  three  terms  of  (106)  are  grouped  together  as  a  quantity  f^ 
Specifically, 


Nx  Nx 


Nx  nDn 


fk  *  ,5,  «  *1*  ‘"x’j  (Vi  +  £  V  ‘Vi  <V« 


(107) 


+  E  E  fiU(PkW- 

j-1  *-l  3  *  3  k 
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Substituting  (31)  into  (107)  yields 


E 

E  f 
i-1 


ik 


(108) 


where 


N  N 
x  x 


“x  »DN 


N  N_„ 

DN  DN 

+  E  Z  “i^Wl 

j-1  «,-!  j  3  k 


(109) 


The  term  multiplying  the  frequency  is  denoted  as  —  . 

\ 

"d»  .  “e  1 

Jl  (<Vl  '  IJ,  «U  "i  »<  11  • 


Specifically, 


(110) 


Although  the  denominator  of  (110)  ia  not  exactly  in  the  form  of  (89), 
it  will  be  shown  that  it  causes  no  problems  in  the  formulation. 

The  quantity  m^  is  certainly  non-negative.  Suppose  the  constraint 
(88)  is  in  a  "less  than"  form,  multiplying  both  sides  by  m^  after 
squaring  (88)  yields 


w]  -  \(u$)2 


(111) 


A  negative  Inverse  operation  preserves  the  inequality,  thus 


-  1 


<  -  1 


\“k. 


\(“?)2  ’ 


(112) 


Adding  f^  to  both  sides  of  (112)  and  employing  (106)  yields 


0  <  f. 


(113) 


\Hy 
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where  N  refers  to  the  number  of  load  cases.  For  simplicity,  the 

L 

formulation  will  continue  with  ■  1,  thus  avoiding  double  subscripting. 
From  (117)  ,  Che  Kuhn-Tucker  optimality  conditions  are  derived.  (It  is 
to  be  remembered  that  L  is  a  function  of  (x)  as  well.)  Application  of 
these  conditions  yields 


3L 

9A, 


-  0 


i-1,  2,  . . . ,  N 


(118a) 


9L  n 


j  -  1.  2 . N 


(118b) 


If  y*>0,  then  g*  -  0  or 
If  y*  -  0,  then  g*£0 


1-1.  2 . \ 


(118c) 


If  y*>0,  then  g*  »  0 
If  y*  -  0,  then  g*£0 


or 


i-1,  2 . N_ 


(118d) 


If  y*£0,  then  g^-0 
If  y^«0,  then  g^lo 

If  y^-0,  then  g^-0 
!f  y^-0,  then  g^O 


or 


or 


i— 1 ,  2,  ...,  N. 


DC 


1-1.  2 . »w 


(118e) 


(118f ) 


gl-0 

•x 


I— 1,  2,  ...,  N. 


(118g) 


The  left-hand  side  of  equations  (118a,  b)  involve  first  derivatives  of 
the  constraint  equations.  The  dependence  In  the  element  sites  for  the 
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constraint  equations  is  that  of  (89)  ,  except  for  the  m^  term  in  the 

frequency  constraints,  and  thus  the  first  derivatives  are  fairly  straight¬ 
forward.  For  the  frequency  constraint. 


NDN 

T  - - - 1 - 

lm  i  \ 

[  Z  eHjm.A.]2 
j-1  J  J 


+ 


(119) 


where  use  has  been  made  of  (108)  and  (110).  Note  that 

gk  +  E  A  \  -  0  (1201 

“  i**l  1  SAj 
thus  preserving  the  form  (94) . 

The  first  derivatives  of  the  constraints  with  respect  to  the 
redundants  must  be  calculated  for  use  in  (118b),  Size  and  frequency 
constraints  are  independent  of  the  loading,  and  thus  independent  of 
the  set  of  redundants.  (For  the  frequency  constraints,  it  is  important 
to  note  the  difference  between  (x),  the  set  of  redundants  as  calculated 
from  the  static  loads,  and  X^,  the  "redundant  mode"  of  the  structure 
when  it  is  in  a  state  of  free  vibration.)  The  first  derivatives  of 
the  other  constraints  are  calculated  as  follows: 

a)  Stress  Constraints  -  With  g  ,  (101)  ,  expressed  as  a 

i 

function  of  the  stress  parameters,  the  derivative  with  respect  to  X^  is 
calculated  by  the  chain  rule  as 
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(121) 


1 


A  O* 
1  i 


3Y. 


9S 


£ 

k-1 


3si.k  s 


where  N  Is  Che  number  of  stress  parameters  for  the  itb  element.  The 
pi 

form  of  ^  is  determined  for  each  element.  (See  (96) 

(99)0.  The  second  set  of  terms  in  the  summation  is  calculated  from 
(29).  Thus, 


3Si.k  .  . i 

bl,kj 


(122) 


b)  Displacement  Constraints  -  Noting  (103)  and  (104), 


3X. 


E  Xk 

JL  2  _ii  . 

A*  k-1  Ak 


(123) 


c)  Equilibrium  Conatrainta  -  Noting  (115)  and  (116), 


3g 

> 

3x, 


Nf  Sk 

ee  hi 

k-1  A, 


(124) 


The  optimization  problem  is  now  specified  explicitly  in  terms 
of  the  unknown  element  sizes  and  redundants,  and  quantities  determined 
from  the  loading,  discretization,  and  choice  of  material  properties 
of  the  elements. 


2.3.3  Description  of  the  Algorithm 

Much  of  the  algorithm  to  be  described  herein  is  a  modification 
of  the  one  developed  in  Reference  4  the  initial  force  method  optimi¬ 
zation  study  program.  «. 


The  guiding  principle  to  be  applied  is  expressed  in  (94)  and  that  equation's 
implicat ions. 

Suppose  a  set  of  element  sizes  are  chosen  such  that  the  set  of 
constraints  (118  c-g)  are  satisfied.  Now,  such  a  set  may  be  quite 
difficult  to  find.  With  (x)  determined  by  >(118g),,  top  priority  may  then 
be  given  to  (118c)  and  (118d) .  The  values  for  the  sizes  and  the 
redundants  are  used  in  (118a)  and  (118b)  to  generate  linear  equations 
in  the  Lagrange  multipliers  u.  With  the  sum  of  the  Lagrange  multipliers 
of  the  type  in  (89a)  as  the  objective  function,  and  >(H8  a,b)  as  the 
constraints,  a  linear  programming  problem  for  the  Lagrange  multipliers  is 
formulated.  The  linear  programming  process  results  in  a  full  vertex  solution, 
i.e.,  for  the  +  N^  equality  constraints  (118  a,b)  1  there  will  be  pre¬ 
cisely  N  +  N  non-zero  valued  y*s  in  the  solution.  It  is  anticipated  that 
L  A 

the  N  non-cons trained  (in  sign)  y  *s  will  be  in  this  basis,  thus  leaving  a 

A  X 

total  of  N  from  among  the  y's,  y  's,  y.'s  and  y  's.  The  non-zero  y’s 
K  A  O  A  U) 

correspond  to  particular  constraints  being  "on",  i.e.,  satisfied  at  equality. 
These  constraints  form  a  system  of  N  +  N  equations  for  the  N  sizes  and 

redundants  for  a  new  design  as  specified  by  the  set  of  non-zero  multipliers. 

This  linear  programming  phase  may  be  initialized  in  one  of  several 
ways.  Firstly,  initial  sizes  may  be  read  from  data  cards;  secondly,  they 
may  be  set  to  minimum  sizes  A^*  by  default;  thirdly,  the  stress  ratio 
method  may  be  employed  either  a  finite  number  of  Iterations,  or  until 
convergence  is  reached.  ThiB  solution  is  known  as  a  fully-stressed  design 
(FSD)  and  is  mathematically  a  full  vertex  solution  when  only  size  and 
stress  constraints  are  considered. 

The  new  design  as  indicated  by  the  set  of  non-zero  y's  is  solved 
from  a  system  of  nonlinear  equations.  If  some  of  the  equations  are  minimum 
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size  constraint  equations,  a  reduction  of  the  system  may  be  accomplished 
by  setting  the  corresponding  sizes  to  the  minimum  directly,  and  eliminating 
these  constraints.  The  remaining  system  is  satisfied  by  employing  a 
Newton-Raphson  procedure.  If  this  procedure  should  not  converge  after  a 
specified  number  of  iterations,  the  results  from  the  last  iteration  are 
used  in  (118a,  b)  and  the  linear  programming  process  is  repeated.  The 
entire  cycle  is  repeated  a  maximum  of  three  times. 

If  the  linear  programming  process  does  not  converge  after  three 
iterations,  an  FSD  design  is  generated.  This  design,  or  the  one  generated 
by  linear  programming  if  that  procedure  converged,  is  used  in  (118a,  b) 
to  determine  values  for  the  corresponding  y'a  to  check  for  their  non- 
negativity,  and  to  check  the  remaining  constraints  (118c-f)  for  violations. 
If  all  the  y*s  are  positive  and  no  constraint  is  violated,  the  solution  is 
found.  If  a  constraint  is  violated,  a  standard  "fix-up”  procedure  is  per¬ 
formed  which  is  dependent  upon  the  particular  constraint  in  question.  For 
instance,  suppose  a  gA  constraint  is  "on"  for  a  given  element  and  the  result¬ 
ing  size  and  sec  of  redundance  yields  a  stress  in  the  element  violating  its 
gQ  constraint.  The  size  is  raised  adequately  to  turn  the  gp  "on"  and  the 
gA  "off".  The  reverse  situation  is  handled  in  a  similar  manner.  If  a 
constraint  is  violated,  all  the  areas  are  multiplied  by  a  common  factor  in 
order  to  satisfy  it.  This  operation,  however,  turns  "off"  the  g^  and  gQ 
constraints.  An  algorithm  to  "fix  up"  g^  constraints  can  be  derived  by 
considering  inverse  mass  to  flexibility  relationships,  to  find  which  elements 
can  have  their  sizes  increased  to  Increase  the  frequency  without  an  excessive 
weight  penalty. 

If  a  particular  y  is  negative,  then  the  corresponding  constraints  is 
turned  "off"  and  the  multiplier  is  set  to  zero. 


S6 


The  current  design  at  this  point  has  no  more  than  N  +  N  constraint 

E  X 

equations  "on".  These  equations,  plus  (118a,  b)  form  a  system  for  the 
element  sizes,  the  redundants,  and  the  non-zero  Lagrange  multipliers 
associated  with  the  constraints.  A  Newton-Raphson  procedure  is  used  to 
solve  this  system.  This  procedure  will  require  the  calculation  of  second 
derivatives  of  the  constraints  which  are  expressed  explicitly  as  follows: 
a)  Minimum  Size  Constraints:  From  (95), 


32gi  =  A* 

2  If  if  1“j=k; 


=  0  otherwise. 


(125) 


The  derivatives  with  respect  to  the  redundants  are  identically  zero. 

b)  Maximum  Stress  Constraints:  From  (101)  it  is  seen  that 


a2  1 

sy\ 


j _ \ 

4  2  V 


if  i=j=k;  *  0  otherwise. 


(128) 


The  other  derivatives  are  calculated  from  (121)  and  (122).  Thus, 


32gi 

°o 

3Xj  \ 


*t  9Y1  .  1 


AiV 


.  ,  3S.  . 

k**l  i,k 


b^  ,  kj  if  i»l;  =  0  otherwise 


(127) 
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A  0  * 
i  i  k-1  m-1 


N_  Np  32 

”i  *i  ° 
l  -  L  1 - - 


as,  ,"as,“  bi  ,kj  bi  ,mi 


i,k  i.m 


(128) 


c)  Displacement  Constraints  -  From  (103), 


a2  * 

8Aj9Ak 


2  -L± 

y  V 


if  j«k;  ■  0  otherwise 


(129) 
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with  is  defined  in  (104);  from  (123), 


3V 


3X 


A 


(130) 


and  the  second  derivatives  with  respect  to  the  redundants  are  zero, 
d)  Frequency  Constraints  -  From  (119), 


.2  < 
3  H 

3Aj3Ak 


(w£)‘ 


ndn 

£ 

1-1 


(IV  2f1i 

\  i  ”  A  3 

[  E  ^mVm  3  J 

m-1 


(131) 


where  the  second  term  is  included  only  if  j«k.  Since  is  independent 
of  the  loading,  all  derivatives  with  respect  to  the  redundants  are  zero. 


e)  Equilibrium  Constraints  -  From  (115), 


<\4.  X  a  _ 

>■  1i  if  j**k;  ■  0  otherwise  (132) 

s*3\  y~ 


where  D 


Ji 


is  defined  in  (116);  from  (124). 


2  i 
g 


(133) 


The  second  derivatives  with  respect  to  redundants  are  zero. 

Assuming  a  given  set  of  equations  form  a  system  that  will  converge 
in  a  Newton-Raphson  procedure,  the  new  design  is  checked  for  feasibility. 

The  violated  constraints  are  "fixed-up"  in  the  manner  prescribed  previously, 
and  those  constraints  with  corresponding  negative  multipliers  are  turned 
off.  The  procedure  of  adding  and  subtracting  constraints  from  the  system 
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to  form  a  new  set  of  equations  for  the  Newton-Raphson  procedure  leads 
eventually  to  a  solution. 

2.4  Rapid  Reanalysis 

The  purpose  of  a  rapid  reanalysis  procedure  is  to  analyze  a  damaged 
structure  using,  as  much  as  possible,  quantities  calculated  in  the  analysis 
of  the  original  structure.  Various  means  to  accomplish  this  task  have 
been  investigated  by  numerous  researchers  References  15  to  20.  .  This 
work,  for  the  most  part,  has  been  based  on  the  matrix  displacement  method 
and  iterative  schemes.  Past  investigators  have  defined  damage  models  by 
removing  structural  finite  elements  entirely  or  reducing  values  of  the 
design  parameters  which  effect  the  structures  flexibility  and  mass. 

J.  S.  Arora,  in  Reference  15  for  example,  defines  the  damage  condition 
as  follows:  "A  damage  condition  for  the  structure  is  defined  to  consist 
of  complete  or  partial  removal  of  selected  members  or  parts  of  a  structure. 
Some  joints  of  the  structure  may  be  removed  as  a  result  of  damage."  These 
authors  applied  this  definition  to  a  helicopter  boom  design  made  of  truss 
members  which  was  subsequently  damaged  due  to  munition  blast  loads  occurring 
inside  or  near  the  boom.  Another  study  illustrates  a  different  approach 
to  the  damage  question.  D.  S.  Scott,  et.al.,  Refere  ce  16  ,  were  concerned 
with  lifting  surface  drag  due  to  holes  in  the  surfaces  caused  by  ballistic 
penetrations.  Their  main  concern  was  the  derivation  of  aerodynamic  loads 
on  surfaces  containing  holes  and  the  consequent  effect  on  aeroelastic 
behavior . 

The  work  of  F.  G.  Hemming,  et.  al..  References  17  and  18, 


considered  battle  damage  or  initial  flaw  propagation  by  removal  of  entire 
finite  elements  or  by  reducing  the  properties  of  the  damaged  elements. 


Property  reduction  could  correspond  to  loss  of  strength  due  to  crack 
propagation  for  instance.  In  either  case  only  a  few  finite  elements  were 
used  to  describe  the  damage  condition.  Additional  work  by  V.  B.  Venkayya,  Re 
erences  19  &  20 j  considered  "severe  damage"  to  a  three  spar  -  five  rib  delta 
wing  by  removing  two  top  and  two  bottom  membrane  elements  and  one  shear 
panel  in  the  midspar.  The  work  being  accomplished  by  the  University  of 
Dayton  Research  Institute,  Dayton,  Ohio,  Reference  21  ,  also  encom¬ 
passes  the  type  of  damage  models  described.  It  was  concluded  from  these 
and  Bell’s  studies  that  for  the  contractual  requirements  two  damage  models 
can  be  defined;  namely,  Types  A  and  B.  Type  A  damage  consists  of  complete 
removal  of  finite  elements  due  to  ballistic  damage  for  instance.  Type  B 
damage  considers  the  reduction  of  finite  element  properties  such  that 
flexibility  is  increased.  This  type  of  damage  could  represent  Increased 
flexibility  due  to  fatigue  cracks  or  small  holes  caused  by  ballistic 
impacts . 

The  rapid  reanalysis  method  developed  for  these  two  types  of  damage 
is  described  below.  A  few  introductory  coimnents  follow.  Note  that 
previously  defined  force  method  matrices  are  used  throughout  the  development . 

Damage  will  be  measured  on  an  element  by  element  basis,  and  will 
consist  of  two  separate  measurements:  d^*,  the  "stiffness  damage"  to  the 

1  and  d^  lie 

between  0  and  1  inclusive,  and  represent  a  fractional  decrease  in  load 
carrying  capacity  and  mass,  respectively.  Generally,  d^  -  0  except  in 
cases  of  physical  removal  of  the  element,  in  which  case  d^1  «  1.  It  is 
conceivable  that  in  cases  of  phase  change  at  high  temperatures  that  d^1 
may  lie  strictly  between  0  and  1,  but,  as  will  be  seen,  the  value  of 


itl1  element,  and  d^*, 


the  "mass  damage".  The  numbers  dR 
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will  not  effect  the  method  of  analysis.  The  value  of  dK*  may  have  any 
value  in  the  range. 

The  following  section  describes  the  procedure  to  be  used  when 
d^*  <  1.  When  dK*  =  1,  a  special  method  must  be  used,  and  it  is  described 
in  Section  2.4.2.  Finally,  the  residual  elastic  strength  of  the  structure 
is  defined  and  discussed. 

2.4.1  Small  Scale  Damage 

Small  scale  damage  is  defined  as  less-than-total  stiffness  damage, 
i.e.,  when  dR*  <  1.  The  effect  upon  element  flexibility  is  thus 

tf±3d  =  [f^/d  -  d^)  (134) 

where  the  d  refers  to  values  in  the  damaged  state.  Note  that  the 
flexibility  increases  with  increasing  damage.  The  difference  between  the 
new  and  old  states  is 

CAfil  =  [fi]d  -  [^3  (135) 

for  each  of  the  Np  damaged  elements. 

Since  the  geometric  layout  of  the  elements  and,  for  the  statics 
case,  the  loading  is  assumed  not  to  change,  the  values  of  [b^],  [d]  and 
{p}  remain  the  same.  Thus,  if  the  definitions 

[<j>3d  -  f>]+[A<J>];l>]d  -  W+[^];[n3d  -  [n]+M  (136) 

are  made,  then 

NE  ne 

[A*]  =  I  [b  1]T[AfiKb11]j[Art  *  E  [b/fCAf.HD.LM 

i-1  i-1  1 


l  [D.3  [Af  ][0  ]. 
i-1 


(137) 
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In  Che  actual  numerical  procedure  che  summation  is  carried  out  only  over 
the  damaged  elements. 

The  equations  of  statics  for  the  damaged  case  are  those  of  (2.20), 
where  d  subscripts  would  be  added  to  all  quantities  used  in  those  equa¬ 
tions,  except  for  (P}.  If  it  is  noted  that 

(X).  -  {x}  +  (AX>  ;  {A}'  -  {A}  +  {aA}  (138) 

a  a 

then,  after  subtracting  (20a)  for  the  undamaged  case  from  (20a)  of 
the  damaged  case,  {AX}  may  be  obtained  as 

{AX}  -  -M”1  (CA<J>1(X}  +  [Ai|»]{P}).  (139) 

The  inverse  of  may  pose  a  problem.  Noting  (136)  and  formulating 

the  expansion 

u]"1  *  ur1  -  urtAoM"1  ♦  tuo) 

will  avoid  the  costly  inverse  operation,  since  [(f)]-1  is  probably  obtained 
in  the  basic  analysis.  The  rate  of  convergence  of  (140)  will  likely 
depend  on  individual  values,  as  well  as  the  proportion  of  elements 
affected  by  damage.  In  order  to  assure  convergence  for  values  of 
above  some  cutoff,  say,  .4,  a  scheme  is  derived  whereby  the  damage  is 
"compounded"  at  a  particular  rate  for  a  finite  number  of  times  until  the 
total,  desired  damage  level  is  achieved.  If  c^*  is  this  compounding  rate, 
then 

1 

ck4  -  l-U-d^)”  (141) 
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where  N  is  an  integral  number  of  steps.  A  criteria  for  choosing  N  would 
be  the  minimum  number  of  compounding  steps  to  achieve  d^™3*,  the  maximum 
damage  level  from  among  all  the  damaged  elements,  with  a  compounding  rate 
of  less  than  .4.  Manipulating  v  (141)  with  d^1  «  .4  yields, 

N  -  log  (l-dkmax)  /  log  .6  .  (142) 

The  value  N  is  then  raised  to  the  nearest  integer  and  used  in  (141)  to 
determine  individual  compounding  rates. 

Let  j-1  iterations  in  the  scheme  be  complete.  The  jth  iteration 
consists  of  deriving  the  new  flexibilities 

tfi]dd)  -  /  d-c/)  (143) 


from  which 


[Af1]<j)  =  [fi]<j)  -  [f1]^'1>  (144) 

is  defined.  As  in  (137)  ,  [A$]  ^  is  defined  as 

NE 

-  Z  [b1l]T[Af.]<J)[bti]  (145) 

i-1  1  1  1 


and  thus 


(146) 


Using  (146)  in  (140)  yields 

(U]JJ'1))"1  -  U*](J)  (a*]*3'1*)'1  +  ...  W) 
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Equations  (143)  -  (147)  are  repeated  for  j-1,  2,  . ...  N.  Note  that 

the  zero**1  state  Is  the  undamaged  state  and  the  NC^  state  is  the  total  damaged 

state. 

Returning  to  the  basic  re-analysis  procedure,  equations  (20b)  for 
the  undamaged  state  are  subtracted  from  the  corresponding  equations  of  the 
damaged  state  to  yield 

([t|»]  +  [At|»])T{AX}  +  [A*]T(X)  +  [AG]{P}  =  (aA)  (148) 

thus  solving  the  re-analysis  problem  completely.  Another  approach  is  to 
find 

Md  '  Md  -  <149> 

and  obtain 

[*]d{P}  -  (A}d.  (150) 

The  global  flexibility  matrix  is  required  when  a  dynamic  re-analysis  is 
needed.  For  this  case,  the  mass  damage  must  also  be  specified.  Its  effect 
is 


<"!>„  *  "i  <151> 

thus 

An^  -  -  dml  m±.  (152) 

Using  (152)  in  (48)  yields 
NE 

i-i  V*!  <153> 
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and  thus 


[M£]d  -  [Mj]  +  L^]  (154) 

where  [AM^]  is  the  proper  collection  of  terms  (153) .  The  matrix 

[QL  is  now  formed  as  in  v  (45)  with  use  of  (154)  and  (149)  and 
a 

eigenvalues  and  eigenvectors  are  now  calculated.  If  they  are  found  using 
some  iterative  technique,  such  as  the  power  method,  the  eigenvectors  of 
the  original  analysis  are  used  as  initial  guesses. 

2.4.2  Large  Scale  Damage 

Large  scale  damage  is  defined  as  the  state  when  at  least  one  -  1. 
This  would  lead  to  an  infinite  flexibility  according  to  (134).  The 
resulting  computational  hazards  would  then  make  the  entire  small  scale 
damage  algorithm  useless.  Note  that  this  is  unlike  the  displacement  method 
whereby  the  damaged  element's  stiffness  is  simply  "removed"  from  the  global 
system. 

An  approach  is  developed  whereby  "total  loss  of  load  carrying  capacity" 
is  translated  into  "carrying  no  load".  That  is,  for  the  ith  damaged  element 

{Si}d  *  {0}  *  (155) 

For  the  collection  of  the  totally  damaged  elements, 

{Sd>d  "  0  (156) 

where  the  d  within  the  brackets  denotes  those  stress  parameters  which  are 
to  be  constrained  to  zero.  Noting  that 

<Sd}  -  [b1d](x}  +  [Dd]{P}  (157) 
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equations  (1)  can  be  extended  by  use  of  Lagrange  multipliers  (X)  as 

6U*  +  6V*  +  6({X}T{S  } .)  -  0.  .  (158) 

a  a 

It  should  be  noted  that  not  only  is  the  geometric  layout  the  same,  and 
thus  [b^]  and  [d],  but  that  the  flexibilities  of  Che  damaged  elements  are 
not  affected,  and  thus  [$],  [\Jj]  and  [0]  remain  the  same.  Taking  varia¬ 
tions  with  respect  to  (x)  and  (P)  in  (158)  yields 


[*]  {X>d  +  [*]{p>  +  [b1d]T  {X}  -  {0} 

(159a) 

W1  (x)d  +  Cfi]{p)  +  [Dd]T{X}  -  {A}. 

(159b) 

Variations  with  respect  to  the  multipliers  return  the  constraints  (156) 
expressed  as 


[b^]  (X}d  +  [Dd]{P)  -  (0).  (160) 

For  the  statlca  case,  the  solution  method  for  (159)  -  (160)  proceeds 
by  introducing  {AX}  and  {AA}  as  in  the  small  scale  damage  case.  Subtract¬ 
ing  (20a)  from  (159a)  yields 

[♦]  (AX)  +  [b^]1  (X)  -  {0}  .  (161) 

The  vector  (AX)  is  solved  in  terms  of  (X),  and  used  in  (160)  to  obtain 

Cb^Kx)  -  [b^iurtb^ftx)  +  [Dd](P)  -  (0).  (162) 

Noting  (157)  ,  the  multipliers  are 

{x}  -  (Cbld3C4]"1Cb1d3T>“1  (sd>.  (163) 
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Subtracting  (20b)  froa  (159b)  yields 

[lj»]T{AX}  +  [Dd]T{X)  -  (A A)  (164) 

An  alternative  approach  is  to  find  CJ^]*  This  will  lend  itself  readily 
to  adoption  for  dynamics.  Solve  (159a)  for  (Xd>,  substitute  into 
(160)  to  solve  for  {A}  and  back  substitute  into  (159b)  to  give  the 
form 

CJL  tP)  =  (A)  (165) 

a 

where 

»d]  -  m  +  ([Dd]  -  Cb1d][4.]‘1M)T(Cb1d][4>]"1[b1d]T)‘]([Dd]-[b1d]U]'1[^). 

(166) 

The  procedure  for  determining  [M^].  is  the  same  as  before,  as  well  as 
the  procedure  for  determining  new  modes  and  frequencies. 

This  method  for  large  scale  damage  is  straight  forward  in  approach 
and  calculation.  One  difficulty  that  may  exist  involves  the  inverse  of 
the  matrix  product  [b^d]  [$]  ^[b^d]^.  If  the  number  of  parameters  in  (sd  J 
exceeds  the  number  of  redundants,  then  this  matrix  is  singular.  The 
physical  interpretation  can  be  either  the  lack  of  rigid  body  motion 
constraints  (as  in  a  rectangular  truss  frame)  or  by  what  will  be  termed 
"node  removal".  For  instance,  if  a  rectangular  plate  is  divided  into  tri¬ 
angular  elements,  and  the  corner  of  the  plate,  represented  by  one  triangle, 
is  removed,  the  corner  node  is  no  longer  physically  there.  The  two 
equilibrium  equations  at  that  node  are  identically  satisfied  (assuming  no 
loads  originally  applied  there)  and  should  be  removed  from  the  analysis. 

A  generalized  method  to  overcome  this  and  similar  difficulties  is  now 
presented. 
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Equations  (14)  are  partitioned  as 


{p}  -  £\JAd]  f  su' 

_'SV 


(167) 


where  the  u  subscript  refers  to  the  undamaged  portion  of  the  structure. 

Perform  a  Gauss-Jordan  elimination  scheme  on  Fa  1  with  column  as  well  as 

*■  uJ 

row  pivoting.  It  is  quite  possible  (and  definite,  if  the  number  of  rows 
of  CAu]  exceeds  the  number  of  columns)  that  this  procedure  yields  the  form 


uu 


»  Aud 


0  |Add 


(168) 


where  the  prime  indicates  that  the  values  are  transformed  from  the  original 
values.  The  lower  partition  of  (168)  are  relations  strictly  among  the 
{s^} .  Note  further  that  unless  {Pj**}  *  (0),  a  contradiction  exists  and 
the  re-analysis  may  not  continue.  Physically,  this  would  imply  that  a 
non-damaged  truss  member,  say,  would  be  required  to  support  a  transverse 
applied  load.  This  would  cause  a  basic  remodeling  of  thie  damaged  structure 
for  analysis,  and  thus  a  rapid  re-analysis  cannot  be  done.  Thus,  for  a 
properly  posed  large  scale  damage  problem,  the  bottom  partition  is 


fs„>  -  (0)  . 


(169) 


A  Gauss-Jordan  elimination  procedure  with  column  pivoting  is  performed  on 
[Add'],  which  should  have  at  least  as  many  columns  as  rows.  This  yields 
the  form 


Cl  !  Add] 


*5 

r»: 


{0} 


(170) 
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The  form  of  (170)  reveals  chat  if  (S  *}  ■  {0i,  then  (Sj)  *  (0),  and 

a  a 

thus  ( S , }  »  {0}  simply  by  equilibrium.  Thus,  in  general,  only  the 
a 

elements  of  {S  *}  need  be  explicitly  set  to  zero  as  in  (156)  in  order 

Q 

to  obtain  the  results  for  the  desired  damage  case.  Furthermore,  if 
every  element  was  removed,  then  [a^"]  “  [a],  and  thus  [a^*]  has  a 


number  of  column  equal  to  the  number  of  redundants,  and  thus  [b^][<j>]  ^[b^]^ 
is  the  same  size  as  [$]. 

Finally,  it  should  be  noted  that  if  some  elements  have  dk*<  1  and 
others  have  d^*  »  1,  that  a  small  scale  analysis  can  be  done  first  with 
the  less  than  totally  damaged  elements,  and  then  that  solution  be  used  as 
the  undamaged  structure  for  the  large  scale  damage. 


2. A. 3  Residual  Elastic  Strength 

The  previous  discussion  has  centered  on  the  results  of  methods  derived 
for  the  rapid  reanalysis  of  damaged  structures.  Output  of  these  methods 
are  such  that  the  vulnerability  of  the  structure  can  be  determined,  this 
being  defined  as  the  ability  of  the  structure  to  withstand  damage  as  mea¬ 
sured  by  residual  elastic  strength.  This  parameter  will  identify  critical 
damage  areas  and  attendant  structural  loads  which  can  cause  complete  failure 
of  the  structure.  The  damage  in  a  structure  may,  for  example,  be  in  the 
form  of  a  crack  in  a  stringer  extending  into  an  adjacent  sheet.  For  analysis 
purposes  the  damaged  portion  of  the  structure  (such  as  the  cracked  plates 
and  portions  of  the  cracked  stringer  in  our  example)  will  be  removed  from 
the  finite  element  idealization  and  an  elastic  reanalysis  performed.  In 
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general,  If  the  structure  had  been  optiaized  and  the  yield  stress  criterion 
was  invoiced  in  the  design,  then  obviously  the  structure  is  weaker  when 
damaged  and  for  the  applied  loads  for  which  the  structure  was  optiaized 
there  will  be  no  residual  strength  up  to  the  yield  stress  values.  However, 
for  loads  less  than  those  used  in  the  optimization  process  there  may  be 
some  residual  strength  remaining  after  damage  which  is  proportional  to  the 
difference  between  the  allowable  yield  stress  and  the  applied  stresses. 

The  strength  of  the  damaged  structure  will  depend  on  the  type  of 
damage  as  related  to  the  failure  mode.  Thus,  for  statically  determinant 
structures  damage  by  removal  of  a  structural  element  can  result  in  a  catas¬ 
trophic  failure.  However,  for  a  redundant  structure  a  damaged  structural 
element  will  cause  a  load  shift  to  other  elements  and  failure  will  be 
governed  by  the  strength  of  remaining  elements.  It  should  be  noted  that  in 
a  highly  redundant  structure  maximum  residual  strength  may  depend  on  the 
failure  of  more  than  one  structural  element. 

In  the  light  of  the  above  complexities  and  in  an  effort  to  estimate 
the  remaining  strength  of  a  damaged  structure,  residual  elastic  strength 
will  be  taken  to  occur  at  loads  which  produce  gross  elastic  stresses  which 
are  equal  to: 

a)  The  yield  stress  of  the  material 

b)  Allowable  stress  used  in  the  optimization  cycle 

c)  The  critical  gross  (nominal)  catastrophic  fracture  stress  based 
on  fracture  mechanics  concepts. 

As  part  of  the  optimization  process,  allowable  stresses  of  the  material 
are  used  as  failure  criteria  in  the  optimum  design  of  the  structure.  Using 
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these  already  identified  allowable  stresses,  the  residual  strength  of  the 
structure  corresponding  to  modes  of  failure  can  be  obtained. 

2. A. A  Computer  Code  &  Illustrative  Examples 

The  rapid  reanalysis  methodology  described  above  was  programmed  into 
a  remote  terminal  interactive,  (VSPC)  computer  code  which  is  aptly  described  and 
illustrated  in  Appendix  B  of  Volume  II  -  User's  Manual.  The  reader  is  referred 
to  this  volume  for  particular  program  details. 

The  use  of  this  computer  program  is  illustrated  here  to  assess  the 
damage  and  residual  elastic  strength  of  three  and  ten  bar  truss  structures. 
Program  input  consists  of  the  characteristics  of  the  optimized  structure  and 
its  loading.  The  damage  level  is  then  expressed  in  terms  of  the  previously 
defined  parameters  d*  and  d*. 

t 

Results  for  a  three  bar  truss  are  shown  on  Figure  7  for  four  damage 
levels.  (Note  that  element  one  sustains  the  damage.)  The  tabular  data 
shows  the  increased  stress  levels  each  bar  experiences  as  damage  levels  in¬ 
crease.  It  Is  observed  that  the  first  element  will  yield  under  the  loads 
used  in  the  optimization  cycle  as  the  damage  level  nears  20%.  This  denotes 
the  vulnerability  of  this  particular  structure  to  the  imposed  damage 
conditions. 

Figure  8  depicts  a  ten  bar  truss  which  was  examined  under  a  variety 
of  damage  cases:  Types  A  and  B  damage  were  evaluated  in  this  instance. 

This  truss  has  been  the  subject  of  study  by  many  investigators  of  weight 
optimization  methods  but  vulnerability  studies  have  only  been  conducted  by 
Messrs.  Venkaya  and  Khot  in  Reference  20.  Thus,  it  provides  a  good  basis 
of  comparison  between  the  Force  and  Displacement  methods  of  structural  analysis. 
Results  of  analyses  conducted  are  given  in  Table  1.  Each  damage  case  is 


given  and  appropriate  results  for  that  case  are  tabulated.  Removal  of  the 
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Figure  7  Three  Bar  Truss  -  Damage  Assessment 
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OPTIMIZATION  CONSTRAINTS: 


TRUSS  CHARACTERISTICS: 


A  .  »  .101n2, 

•  in  * 

*  -  t  25,000  p«l 


2.0  In, 


p  -  .10  lbs/in3 
E  -  1.0  x  107  ps 


RESULTS  OF  OPTIMIZATION:  -  5062  lb«. 
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.57 

2 

.10 

7 
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22.73 

8 

21.33 

4 

15.45 

9 
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5 

.10 

10 

.  10 

Figure  8  Ten  Bar  Truss 
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TABLE  1  TEN  BAR  TRUSS  -  DAMAGE  ASSESSMENT 


DAMAGE 

CASE 

<*k 

*aax^ 

Aaax 

COMMENT 

1,1-9 

0 

6525 

*5-24970 

-2.0 

Not  Critical 

.5 

13014 

°S-25149 

-2.5 

Not  Critl-al 

.  1 

63040 

*9-63040 

-6.1 

Structure  Collapse 

1.0 

- 

*10>1. 0E6 

-146.4 

Structure  Collapse 

2,1-10 

1.0 

- 

*5-25050 

-1.99 

Not  Critical 

3,1-7 

1.0 

- 

"S-800,000 

1 

w 

o 

o 

Structure  Collapse 

4,1-6 

1.0 

- 

*5-25050 

esi 

Not  Critical 

1  ■  Daaaged  Bar  Element  Number 

a  In  pounds  par  aquare  Inch 

A  in  lnchaa 

Aaax  *  Displacement  of  Structure 
VaaXj-  Max.  Straaa  Occura  In  Jtb  Eleaant 
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V 


ninth  element,  indicated  by  i*9,  produces  a  condition  of  complete  collapse  of  the 
structure  as  evidenced  by  a  maximum  displacement  of  146.4  in.  and  stress  levels 
exceeding  one  million  pounds  per  square  inch.  It  must  be  noted  that  this 
result  is  in  agreement  with  Ref erence  20  and  unlike  the  reanalysis  method 
used  in  that  reference  no  iteration  procedure  is  used.  The  additional 
conditions  tabulated  illustrate  the  vulnerability  of  the  truss  to  damage 
and  demonstrates  the  importance  of  finding  those  members  which  cause  the 
entire  structure  to  fail.  Conditions  such  as  these  permit  definition  of 
residual  strength  levels. 

Figure  9  displays  the  approach  to  determining  the  elastic  residual 
strength  of  the  three  bar  truss  through  use  of  the  method  discussed  above. 

The  stresses  in  each  element  of  the  optimized  truss  were  limited  to  i25,OOQ 
psi.  Assuming  that  a  60%  damage  level  (d^=,6)  will  be  sustained  by  the 
first  element  then  the  maximum  load  that  the  damaged  truss  can  accommodate 
is  Pd».40Po,  where  PQis  the  loading  that  the  optimized  structure  was  designed 
to.  In  this  instance  the  residual  elastic  strength  level  is  40%. 
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3.0  OPTFORCE  II  PROGRAM 


3.1  Efficiency  Studies 

Task  VI  of  the  work  effort  required  that  the  efficiency  of  the  force 
method  optimization  program,  0PTF0RCE  II,  be  compared  with  an  existing  dis¬ 
placement  method  optimization  program  OPTIM  III  (Ref.  10).  Much  of  the 
optimization  literature  in  which  solution  methods  are  compared  report  the 
number  of  analyses  required,  apparently  using  this  as  a  measure  of  effi¬ 
ciency  of  the  method  being  discussed.  This  may  be  a  valid  measure  of 
efficiency,  but  it  appears  to  be  dependent  on  specific  machine  capability. 
For  example,  if  only  a  few  analyses  are  required  but  each  analysis  is 
rather  complex,  requiring  somewhat  more  computer  time,  this  method  would 
appear  to  be  more  efficient  on  the  basis  of  number  of  analyses  required 
than  a  method  requiring  many  analyses,  each  of  very  short  duration.  On 
the  other  hand,  using  computer  time  (seconds)  as  a  measure  of  efficiency 
requires  that  the  problems  to  be  used  for  comparison  be  exercised  on  the 
same  computer  system  to  avoid  faster  or  slower  machine  influences.  The 
efficiency  studies  discussed  herein  were  based  on  computer  time  but  due 
note  was  taken  of  the  number  iteration  cycles  as  well.  Each  of  the  problems 
were  solved  either  on  Bell's  IBM  3031  or  IBM  3033  computer  and  the  time 
recorded  is  the  cpu  computer  time,  in  seconds,  to  optimize  the  structure. 

It  became  very  evident,  as  the  study  progressed  and  structural 
optimization  problems  were  solved,  that  different  design  variable  popula¬ 
tions  were  obtained  for  the  same  minimum  weight  value  using  both  of  the 
above  computer  codes  and  other  reference  material.  Thus,  the  "efficiency" 
studies  were  expanded  to  include  program  accuracy  as  well.  The  "accuracy" 
of  a  particular  optimization  code  was  gauged  by  how  well  it  predicts  the 
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crue  minimum  weight  solution.  The  measure  used  to  accomplish  this  task 
was  the  agreement  obtained  between  different  codes  and  known  analytical 
solutions.  The  ensuing  discussion  considers  each  of  the  above  elements  to 
assess  the  efficiency  of  the  OPTFORCE  II  code.  Three  structures  were 
optimized  and  are  discussed  below  in  turn. 

3.1.1  Seventeen  Bar  Truss 

The  configuration  of  the  seventeen  bar  truss  is  given  in  Figure  10 
with  two  external  loading  cases;  each  load  is  applied  simultaneously,  they 
are  not  considered  as  multiple  load  cases.  Material  properties  and  con¬ 
straints  for  each  loading  case  are  given  in  Table  2.  Only  stress  and 
minimum  size  constraints  were  considered.  The  OPTFORCE  II  NASTRAN  compatible 
input  data  for  this  structure  is  displayed  in  Figure  11  for  Case  2.  The 
OPTIM  III  input  data  file  is  essentially  identical  to  those  data  shown  in 
Figure  11. 

Case  1  results  using  OPTFORCE  II  and  OPTIM  III  optimization  programs 
are  displayed  in  Table  3  and  Figure  12.  Data  from  Drs.  Khot  and  Berke's 
research.  Reference  22,  are  also  included  for  comparative  purposes.  Their 
solution  procedures  are  based  on  the  use  of  two  algorithms.  The  first 
algorithm  Is  a  recurrence  relation  based  on  the  fully  stressed  design 
(FSD'  criterion  and  the  second  algorithm  uses  a  recurrence  relationship 
based  on  equivalent  displacement  constraints.  Examination  of  the  results 
shows  the  excellent  agreement  obtained  among  the  methods  used.  Note  that 
OPTIM  III  and  Khot  and  Berke's  methods  are  "displacement"  based  whereas 
OPTFORCE  II  is  "force"  based.  Of  particular  interest  is  the  initial  weight 
value,  ,  and  number  of  iterations,  I,  comparisons.  This  particular 
application  of  OPTFORCE  II  used  the  stress  ratio  option  to  arrive  at  an 
initial  guess  for  the  design  variable  vector.  This  procedure  yielded 
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Grid  Point 

Axis 

Load  i 

9 

3, 5, 7, 9 

+Y 

+Y 

-100000.0  lbs. 
-100000.0  lbs. 

Figure  10  Seventeen  Bar  Truss 
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TABLE  2  MATERIAL  PROPERTIES  &  CONSTRAINTS  -  SEVENTEEN  BAR  TRUSS 


(1)  Material  Properties 
Aluminum: 

E  =  30x10**  psi 

p  *  .268  lbs.  per  cubic  inch 
v  »  .30 

(2)  Minimum  size  (size  constraints) 

2 

Bar  element  area  *  .10  in 

(3)  Allowable  stress  (stress  constraints) 

Case  1: 

0^  *  -50000.0  psi,  0u  -  50000.0  psi  all  elements 
Case  2: 

o.  ■  -50000.0  psi,  a  -  50000.0  psi  all  elements 
except  No.  2,  6,  10  U 

o.  -  -125,000.0  psi,  0  ■  125,000.0  psi  elements 
No.  2,  6,  10  u 


I 


10  TITLE 

11  SEVENTEEN  BAR  TRUSS-FOUR  LOADS 


20  GRID 

1 

0.0 

0.0 

0.0 

123 

30  GRID 

•y 

0.0 

100.0 

0.0 

123 

40  GRID 

3 

100.0 

0.0 

0.0 

3) 

SO  GRID 

4 

100.0 

100.0 

0.0 

3 

60  GRID 

S 

200.0 

o.a 

0.0 

3 

70  GRID 

6 

200.0 

100.0 

0.0 

3 

80  GRID 

7 

300.0 

0.0 

0.0 

3 

90  GRID 

8 

300.0 

100.0 

0.0 

. 

3 

100  GRID 

9 

400.0 

0.0 

0.0 

3 

110  OPLOADS 

1 

i 

150  OPTIH 

YES 

NO 

YES 

OPT 

160  FORCE 

1 

3 

1.0+5 

0.0 

-1.0 

0.0 

170  FORCE 

1 

5 

1.0+3 

0.0 

-1.0 

0*0 

180  FORCE 

1 

7 

1.0+5 

0.0 

-1 .0 

0.0 

190  FORCE 

1 

9 

1.0+5 

0.0 

-1.0 

0.0 

200  CONROD 

1 

2 

4 

1 

0.10 

0 

0 

o 

210  CONROD 

2 

2 

3 

2 

0.10 

0 

0 

o 

220  CONROD 

3 

1 

3 

1 

0.10 

0 

0 

o 

230  CONROD 

4 

3 

4 

1 

0.10 

0 

0 

o 

240  CONROD 

S 

4 

6 

1 

0.10 

0 

0 

o 

2S0  CONROD 

6 

4 

5 

3 

0.10 

0 

0 

0 

260  CONROD 

7 

3 

5 

1 

0.10 

0 

0 

o 

270  CONROD 

8 

S 

6 

1 

0.10 

0 

0 

o 

280  CONROD 

V 

6 

8 

1 

0.10 

0 

0 

o 

290  CONROD 

10 

6 

7 

4 

0.10 

0 

0 

o 

300  CONROD 

11 

3 

7 

1 

0.10 

0 

0 

o 

310  CONROD 

12 

7 

8 

1 

0.10 

0 

0 

o 

320  CONROD 

13 

a 

9 

1 

0.10 

0 

0 

o 

330  CONROD 

14 

7 

9 

1 

0.10 

0 

0 

o 

340  CONROD 

IS 

1 

4 

1 

0.10 

0 

0 

o 

350  CONROD 

16 

3 

6 

1 

0.10 

0 

0 

o 

360  CONROD 

17 

S 

8 

1 

0.10 

0 

0 

0 

370  HAT1 

1 

30. +6 

0.3 

0.268 

+NATA 

380  +HATA 

-5.0+4 

♦5.0+4 

390  NAT1 

2 

30. +6 

0.3 

0.268 

+NAT8 

400  +NATB 

-12.344 

+12.5+4 

410  HAT1 

3 

30. +6 

0.3 

0.268 

+NA7C 

420  +HATC 

-12,5+4 

+12.5+4 

430  MAT1 

4 

30. +6 

0.3 

0.268 

+MATD 

440  +NATD 

-12.S+4 

+12.5+4 

Flugre  11  Opt force  11  Input  Dnta  -  Seventeen  Bar  Trues,  Case  2 


TABLE  3  SEVENTEEN  BAR  TRUSS  RESULTS  -  CASE  1 


i 


-  1315.75  lbs,  a  value  very  near  the  minimum  weight  value  of  W»1295.49  lbs. 
Observe  that  this  is  not  the  case  for  the  other  computer  codes  shown.  The 
OPTIM  III  code  does  not  have  this  capability  and  uses  the  minimum  values  of 
the  design  variables  as  the  initial  guess.  The  OPTFORCE  II  code  yielded 
the  optimum  weight  solution  in  one  iteration  (Linear  Programming  Phase) 
whereas  the  other  methods  required  30  to  36  iterations.  Computer  time 
recorded  shows  that  the  OPTFORCE  II  code  requires  6.55  more  cpu  seconds  then 
the  OPTIM  III  code.  The  inset  on  Figure  12  indicates  that  a  solution  is 
not  feasible  in  the  Linear  Programming  phase.  A  non-feasible  solution 
violates  at  least  one  of  the  constraints.  For  a  given  problem,  it  is  pos¬ 
sible  that  no  feasible  solutions  exist.  Conversely,  a  feasible  solution  is 
any  solution  of  the  constraint  equations. 

Case  2  results  are  shown  in  Table  4  and  Figure  13.  No  agreement 
is  obtained  between  the  OPTFORCE  II  and  OPTIM  III  codes,  however,  agreement 
occurs  between  OPTFORCE  II  and  Reference  22  solutions.  It  is  discussed  in 
that  reference  that  the  FSD  algorithm  yielded  a  "minimum  weight"  design  of 
3081.62  lbs.  which  is  the  same  design  obtained  using  the  OPTIM  III  code.  As 
discussed  in  Reference  22  this  solution  is  a  non-optimum  design  and  is  not 
the  correct  one.  It  was  further  shown  in  Reference  22  that  the  equivalent 
displacement  algorithm  yielded  the  optimum  design,  it  being  of  weight 
equal  to  2460.24  lbs.  Thus  as  seen  the  correct  distribution  of  areas  and 
weight  for  the  optimum  design  were  obtained  by  the  OPTFORCE  II  code  and 
not  the  OPTIM  III  code.  Note  that  the  stress  distribution  for  all  three 
methods  shown  are  identical  which  illustrates  that  there  can  be  more  than 
one  design  with  tl.e  same  stress  distribution  but  different  weights.  This, 
of  course,  was  also  observed  and  quoted  in  Reference  22.  Computer  time 


TABLE  4  SEVENTEEN  BAR  TRUSS  RESULTS  -  CASE  2 


► 


CPU  14.16  sec.  8.78  sec.  Not  given 


Csa'i)  ift 


Figure  13  Seventeen 


displayed  in  Table  4  again  shows  that  the  force  method  incurs  additional 
solution  time  but  requires  one-sixth  the  number  of  iterations. 

3.1.2  Four  Bar  Pyramid 

Figure  14  displays  the  geometrical  configuration  of  the  four  bar 
pyramid.  The  external  loading  applied  to  this  structure  is  also  shown. 

This  loading  system  is  not  considered  as  a  multiple  load  case;  rather  each 
of  the  loads  are  applied  to  the  structure  simultaneously.  Material  proper¬ 
ties  and  constraints  are  listed  in  Table  5.  Only  stress  and  minimum  size 
constraints  were  considered.  OPTFORCE  II  input  data  is  displayed  in 
Figure  15;  OPTIM  III  data  is  essentially  identical  to  that  shown  in  that 
figure.  Analyses  results  are  given  in  Tables  6  and  7  and  Figures  16 
and  17. 

Three  cases  were  considered  differing  only  in  the  value  of  the 
minimum  size  constraints.  The  minimum  weight  for  these  cases  was  calcu¬ 
lated  to  be  65.76  lbs.;  the  value  accepted  as  the  optimum  by  various 
investigators.  Displacement  and  stress  values  obtained  were  identical 
for  all  cases  displayed,  however  distinctly  different  designs  are  evident 
as  shown  by  the  values  of  the  areas  obtained.  Case  1,  OPTFORCE  II  analysis, 
starts  with  an  initial  guess  vector  for  the  design  variables  based  on  the 
stress  ratio  method  which  yielded  the  area  vector  (.43,  1.76,  1.26,  .55). 
This  vector  was  subsequently  revised  through  use  of  the  Linear  Programming 
subroutine  to  that  shown  in  Table  6  and  Figure  16.  Identical  results 
were  obtained  to  those  shown  when  the  stress  ratio  option  was  not  exercised. 
Note,  however,  the  additional  computer  time  and  iterations  needed.  Figure 
16  aptly  portrays  these  results.  OPTIM  III  code  computations  converged  to 

the  area  vector  (.43,  1.76,  1.26,  .55)  a  result  obtained  by  Gellatly, 

2 

(Reference  23),  and  Venkayya  (Reference  24)  using  Amin  ■  0.0  in  .  Reduction 
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TABLE  5  MATERIAL  PROPERTIES  &  CONSTRAINTS  -  FOUR  BAR  PYRAMID 


10  TITLE 


20  FOUR 

BAR 

PRVAMID 

30  GRID 

1 

0.0 

0.0 

0.0 

40  GRID 

2 

204.0 

0.0 

0.0 

50  GRID 

3 

204.0 

192.0 

0.0 

60  GRID 

4 

0.0 

192.0 

0.0 

70  GRID 

5 

60.0 

120.0 

96.0 

80  SPC1 

1 

123 

1 

2 

3 

4 

90  OPTIM 

NO 

NO 

YES 

0 

100  OPLOADS 

1 

1 

110  FORCE 

1 

5 

10000.0 

1.0 

2.0 

143  CONROD 

1 

1 

5 

1 

.20 

0 

150  CONROD 

2 

4 

5 

1 

.20 

0 

160  CONROD 

3 

3 

5 

1 

.20 

0 

170  CONROD 

4 

2 

5 

1 

.20 

0 

180  MAT1 

1 

10.046 

.30 

.10 

190  +HAT  -25000. 025000. 0 
200  ENDDATA 


OPT 


-6.0 

0  0 

0  0 

0  0 

0  O 


MAI 


Figure  15  Optforce  II  Input  Data  -  Four  Bar  Pyramid,  Caae  1 


J 
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TABLE  6  FOUR  BAR  PYRAMID  -  CASES  1  &  2 


Case  1  A  .  ■  . 

min 

20  in2 

El. 

OPT FORCE  II 

OPTIM  III  1 

No. 

Area (in2) 

Stress(psl) 

Area (in2) 

Stress(psi) 

1 

.70 

-25000.0 

.43 

-25000.0 

2 

1.53 

-25000.0 

1.76 

-25000.0 

3 

1.57 

-25000.0 

1.26 

-25000.0 

A 

.20 

-25000.0 

.55 

-25000.0 

w 

i 

65.76/13.94*  lbs. 

122.35  lbs. 

U 

65.76  lbs. 

65.76  lbs. 

m 

1 

X5 

.21  in. 

.21  in. 

Ys 

-.12  in. 

-.12  in. 

Zs 

-.69  in. 

-•69  in. 

CPU 

6.30/23.54*  sec. 

2 

46  sec. 

1 

4/25* 

3 

*  Scress  ratio  option  not  used  for  initial  quess. 


Case  2  A 


min 


.10  in2 


El. 

0PTF0RCE  II 

OPTIM  III 

No. 

Area (in2) 

Stress(psi) 

Area (in2 ) 

Stress(psi) 

1 

.78 

-25000.0 

.43 

-25000.0 

2 

1.47 

-25000.0 

1.75 

-25000.0 

3 

1.66 

-25000.0 

1.26 

-25000.0 

4 

.10 

-25000.0 

.55 

-25000.0 

u 

i 

65.76  lbs. 

122. 

35  lbs. 

W 

m 

65.76  lbs. 

65. 

76  lbs. 

Xs 

.21  in. 

.21  in. 

Ys 

-.12  in. 

-.12  in. 

Z5 

-.69  in. 

-.69  in. 

CPU 

6.71  sec. 

2.40  sec. 

OPT  FORCE  II 

INITIAL  QUESS  -  STRESS  RATIO 
METHOD  USING  MIN.  SIZES,  W-65.76  LBS. 

♦ 

ENTER  LINEAR  PROGRAMING  ITER  NO.  1 
SOLUTION  IS  FEASIBLE.  DESIGN  CONVERGED 
IN  PARTIAL  HRHTON-RAPHSOM 

LAGRANGE  MULTIPLIER  CHECK:  ALL  X'S 
POSITIVE  USE  DESIGN  PROM  LINEAR 
PROGRAMING . 


SOLUTIONS 

•  INITIAL  QUESS  -  MIN.  SIZES  H-1J.94  LBS. 

* 

•  ENTER  LIN  PROGRAMING,  ITER.  NO.  1 
SOLUTION  la  FEASIBLE.  DESIGN  CONVERGED 
IN  PARTIAL  NEVTON-RAPHSON. 

+ 

•  LAGRANGE  MULTIPLIER  CHECK:  NEGATIVE  X  'S 
ENTER  FULL  N-R. 

i 

•  DESIGN  CONVERGED  IN  FULL  N-R. 
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INITIAL 


OPT FORCE  II  SOLUTIONS 


•  INITIAL  QUESS  -  STRESS  RATIO 
METHOD  USING  MIN.  SIZES, 
W-65.76  LBS. 

+ 

•  ENTER  LINEAR  PROGRAMMING, 
ITER  NO.  1.  SOLUTION  IS 
FEASIBLE.  DESIGNED  CONVERGED 
IN  PARTIAL  NEWTON-RAPHSON. 

+ 

a  LAGRANGE  MULTIPLIER  CHECK: 
ALL  X’S  POSITIVE.  USE  DESIGN 
FROM  LINEAR  PROGRAMMING. 


CASE  2.  A 


OPTIM  III 
W  -65.76  LBS 


OPTIM  III 
W-65.76  LBS 


\ 


'ZcD 


‘‘OPTFORCE  II 
W^-65,76  LBS. 


V 


OPTFORCE  II 
W-65.76  LBS. 


i  a  I 


a  a  I 


2  cn 

NO.  ITERATIONS 

S3 

NO.  ITERATIONS 

Figure  17 

Four  Bar  Pyramid 

-  Cases  2  &  3 

2 

of  the  minimum  area  constraint  to  .10  in  gave  the  results  tabulated  for 

Case  2  in  Table  6  and  Figure  17.  In  this  instance  0PTF0RCE  II  computa- 

tions  yielded  the  vector  (.78,  1.47,  1.66,  .10)  and  OPTIM  III  again 

yielded  the  vector  (.43,  1.75,  1.26,  .55).  Further  reduction  of  the  area 

2 

constraint  to  .0001  in  yielded  the  area  vectors  listed  under  Case  2, 

Table  7  and  Figure  17. 

Perusal  of  the  above  results  suggests  that  multiple  sets  of  design 

variables  satisfy  the  minimum  weight  state.  Such  is  the  case  in  this 

instance  as  thoroughly  discussed  by  Gellatly  and  Venkayya.  Their  results 

are  shown  under  Case  4,  Table  7.  Further  work  done  by  Venkayya  showed 

that  multiple  minima  exist  for  the  structure  under  consideration.  He 

states  that  "the  design  weight  of  65.76  lbs.  appears  to  be  the  absolute 

minimum  for  this  structure  and  there  are  three  designs  having  the  minimum 

weight.  The  other  two  may  be  considered  relative  minimums".  The  design 

2 

vectors  listed  in  Reference  24,  for  A  -  0.0  in  ,  are  (0.0,  2.105,  .770, 

min 

1.097)  (.859,  1.406,  1.746,  0.0)  and  (.430,  1.755,  1.258,  .548)  with  the 
latter  associated  with  the  absolute  minimum  weight  listed  in  Table  7 
Comparison  of  these  area  vectors  with  those  tabulated  shows  that  the 
OPTIM  III  solutions,  regardless  of  the  minimum  size  constraint  value, 
agrees  with  the  absolute  minimum  solutions  of  Gellatly  and  Venkayya  using 
Amin  “  The  OPTFORCE  II  solution  procedure  tends  to  one  of  the  relative 

minimums  as  the  minimum  size  constraint  vanishes  which  thoroughly  agrees 
with  the  results  obtained  by  the  above  two  investigators.  Thus ,  it  is 
concluded  that  OPTFORCE  II  yields  the  correct  optimal  solution  for  Case  l; 
the  initial  design  problem.  Note  that  Hsrless,  Reference  25,  obtains 
yet  another  solution.  This  solution  listed  in  Table  7  is  very  close 
to  one  of  the  relative  minimums  obtained  by  Venkayya.  Examination  of 


Tables  6  and  7  shows  that  OPT FORCE  II  requires  more  cpu  seconds  than 
that  required  for  OPTIM  III  for  nearly  the  same  number  of  solution  itera¬ 
tions.  What  is  most  Interesting,  as  was  observed  in  the  seventeen  bar 
truss  solutions,  is  that  the  initial  weight  values  *  V, ,  used  in  the 
OPTFORCE  II  code  are  very  close  to  the  optimum  weight  W  .  This  is  not 

ID 

the  situation  with  the  other  reference  computer  codes.  Thus,  the  initial 
guess  using  the  optional  stress-ratio  method  is  definitely  an  asset  to  the 
User  and  as  demonstrated  in  Case  1  its  use  reduces  computer  cpu  time  and 
the  number  of  iteration  cycles  required. 

3.1.3  Wlngbox 

The  configuration  of  the  eighteen  element  wlngbox  is  given  in 
Figure  18  with  the  location  of  the  single  external  load.  Material  proper¬ 
ties  and  problem  constraints  are  listed  in  Table  8;  only  minimum  size 
and  stress  constraints  were  considered.  Figure  19  displays  OPTFORCE  II 
input  data.  The  OPTIM  III  input  data  file  is  essentially  the  same  at 
that  shown  in  Figure  19. 

Results  of  the  analyses  conducted  are  displayed  in  Table  9  and 
Figure  20.  Examination  of  these  results  shows  that  the  minimum  weights 
obtained  are  relatively  close  to  one  another  with  the  OPTFORCE  II  code 
yielding  a  weight  nearly  10%  less  than  that  calculated  using  the  OPTIM  III 
code.  The  displacement  characteristics  obtained  from  both  codes  were 
essentially  identical.  Viewing  the  design  variable  vector  shows  that 
two  different  designs  were  obtained.  The  most  radical  departure  between 
these  designs  appears  in  the  rod  and  web  finite  elements  aligned  with 
the  applied  force  (elements  4,  5,  6,  10,  11  and  12).  The  quadrilateral 
membrane  elements  also  exhibit  the  same  behavior.  This  disagreement  is 
attributed  to  two  factors;  namely,  1)  difference  in  optimization  solution 
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10  II I IX 

JO  fc  1GHTEEN  ELEMENT  UINGBOX 


30 

OR  IB 

1 

0. 

0 . 0 

1 0 .  3 

•♦0 

GRID 

7 

100/ 

0.0 

1J .  0 

50 

OR  111 

3 

0. 

70.0 

1  0 . 0 

*0 

OR  U< 

4 

100. 

70 . 0 

11.11 

70 

OR  ID 

5 

0. 

140.0 

10. 0 

80 

GRID 

6 

100. 

140.0 

13.0 

90 

GRID 
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Figure  19  Optforce  II  Input  Data  -  Wingbox 
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TABLE  9  WINGBOX  RESULTS 


El .  Element 

No.  Type 


OPTFORCE  II 

Aj*  j  Stress(psi) 

OPTIM  III 

A^*  j  Stress(psi) 

.1000 

-9897.0 

.2404 

-7996.0 

.1000 

-7684.0 

.1010 

-6073.0 

.1000 

-1274.0 

.1010 

-217?. 0 

4.5500 

-10000.0 

.3437 

-8927.0 

2.1060 

-10000.0 

.1446 

-8349.0 

.2527 

-10000.0 

.1010 

-6224.0 

.0200 

.-1 590.0 

.0202 

-1446.0 

.0200 

-731.0 

.0202 

2321.0 

.0200 

5018.0 

.0202 

8052.0 

.0966 

5774.0 

.0449 

9896.0 

.0897 

5774.0 

.0420 

9903.0 

.0865 

5774.0 

.0466 

9918.0 

.0200 

2157.0 

.0202 

5714.0 

.0200 

3429.0 

.0202 

6639.0 

.0200 

4014.0 

.0202 

5153.0 

.0434 

9999.0+ 

.0931 

8735.0+ 

.0326 

10000.0+ 

.0643 

8568.0+ 

.0211 

10000.0+ 

.0303 

8496.0+ 

136. 

64  lbs. 

234. 

69  lbs. 

135. 

96  lbs. 

148. 

15  lbs. 

2.54  in. 

2. 

42  in. 

19.96  sec.** 

.96 

sec.** 

14 

6 

Design  variable  value,  rod  cross-sectional  area  (in  ) 
web  thickness  (in.),  quad,  membrane  thickness  (in.) 

IBM  3033 

Mises-Henchy  stress  criteria  shown. 


Figure  20  Wingbox  Results 


procedure,  e.g.  "force"  method  versus  "displacement"  method  and  2)  possible 
finite  element  formulations.  As  noted  in  the  discussions  of  the  seventeen 
bar  truss  and  four  bar  pyramid  solutions  the  OPTFORCE  II  code  appears  to 
be  the  more  accurate  one  when  the  weight  optimization  option  is  exercised. 
This  fact  is  particularly  emphasized  when  the  statics  option  was  used. 
Results  of  this  exercise  demonstrated  that  both  codes  yielded  nearly  ident¬ 
ical  values  of  the  displacement,  stress  and  reaction  vectors.  Cpu  times 
recorded  for  the  statics  case  were  .88  sec  using  OPTFORCE  II  and  .68  sec. 
for  OPTIM  III  a  negligible  difference.  The  stresses  displayed  in  Table  9 
are  the  element  stresses  themselves  with  exception  of  those  given  for  the 
quadrilateral  membrane  elements.  Three  element  stresses  are  calculated 
(ox,  Oy,  Txy)  for  these  elements  and  these  have  been  combined  using  the 
Mises-Henchy  stress  failure  criteria  for  the  sake  of  brevity  and  comparison. 

Computer  time  recorded  shows  that  OPTFORCE  II  used  19.96  sec.  of  cpu 
time  versus  .96  sec.  for  OPTIM  III  solutions  on  the  IBM  3033  computer.  The 
number  of  iterations  required  for  OPTFORCE  II  are  greater  than  those  shown 
for  OPTIM  III.  Examination  of  Figure  20  and  the  computer  output  shows 
that  the  minimum  weight  value  obtained  from  OPTFORCE  II  has  converged  in 
less  than  the  number  of  iterations  shown  (successive  weight  values  only 
differing  in  the  second  or  third  decimal  place).  Additional  iterations 
were  required  to  obtain  convergence  of  the  design  variable  values.  This 
was  observed  to  be  case  in  other  applications  of  the  subject  computer  codes. 
3.1.4  Concluding  Remarks 

The  above  discussions  clearly  illustrate  the  fact  that  there  can  be 
more  than  one  design  with  the  same  stress  distribution  and  minimum  weight 
but  different  design  variable  values  or  so-called  area  vectors.  It 
appears  that  this  fact,  the  accuracy  component  of  the  efficiency  study. 
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is  the  driving  force  behind  optimization  computer  code  selection  rather 
than  computer  time  and/or  the  number  of  iterations  required.  The  OPTFORCE 
II  code  not  only  yields  the  correct  minimum  structural  weight  but  the 
correct  design  variable  population.  This  conclusion  cannot  be  over¬ 
emphasized  to  the  potential  structural  designer.  Computer  time  is  begin¬ 
ning  to  have  less  relevance  as  a  measure  of  computer  code  efficiency  than 
in  the  past  due  to  the  higher  speeds' attained  with  the  newer  machines. 

Several  of  the  solutions  discussed  in  this  section  were  re-run  using  Bell's 
more  recently  acquired  IBM  3033  computer.  Identical  minimum  weight  solutions 
were  obtained  using  both  OPTFORCE  II  and  OPTIM  III  codes  at  approximately 
one-fifth  the  original  cpu  times  quoted  using  the  IBM  3031  machine.  It  is 
concluded  from  these  facts  and  the  discussions  of  minimum  weight  solutions 
that  the  OPTFORCE  II  code  is  more  "efficient"  than  the  reference  OPTIM  IIT 
code. 

3.2  Applications 

The  development  of  the  structural  optimization  technique  using  the 
force  method,  which  has  been  amply  described  in  Section  2.0,  is  further 
illustrated  below  using  the  swept  wing-box  idealization  shown  in  Figure 
21.  Two  material  cases  were  chosen;  namely,  Case  1  an  all-aluminum 
structure  and  Case  2  a  graphite/epoxy-aluminum  structure.  Vertical  loads 
are  applied  at  grid  points  4,  7,  and  10  and  are  considered  as  a  single 
load  case. 

Case  1  idealization  consists  of  modeling  the  upper  skins  with 
quadrilateral  finite  elements  whereas  the  spars  and  ribs  are  modeled 
using  the  symmetric  shear  panel  elements.  Spar  and  rib  caps  are  idealized 
using  rod  finite  elements.  It  is  noted  that  only  one-half  of  the  wing-box 
structure  needed  to  be  modeled  by  virtue  of  the  geometric  symmetry  of  the 
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ROOT  SECTION  (NO  SCALE) 


Figure  21  Swept  Wingbox  Configuration 


structure  and  use  of  the  symmetric  shear  panel  element.  As  a  result,  only 
six  quadrilateral  membrane,  fifteen  shear  panel  and  fifteen  rod  finite 
elements  are  used  to  idealize  the  entire  wing-box  yielding  a  total  of 
thirty-six  elements.  Twenty-seven  degrees-of-freedom  describe  the  struc¬ 
tural  behavior.  The  design  variables  are  quadrilateral  membrane  thickness 

(t  ),  rod  cross-sectional  area  (A)  and  shear  web  thickness  (t  );  thirty- 
m  w 

six  in  number.  Material  properties;  minimum  sizes  of  design  variables 
and  allowable  stress  levels  are  given  in  Table  10.  NASTRAN  compatible 
input  for  this  structure  is  displayed  in  Figure  22. 

Case  2  idealization  consisted  of  again  modeling  the  upper  skins 
using  quadrilateral  membrane  elements.  However,  in  this  case,  the  skins 
consisted  of  four  layers  of  elements  containing  common  grid  points.  Each 
of  the  layers  contain  a  different  fiber  orientation  in  the  graphite/epoxy 
material.  The  spars,  ribs,  and  caps  were  modeled  as  in  the  Case  1  structure. 
A  total  of  fifty-four  elements  and  twenty-seven  degrees-of-freedom  results 
from  this  idealization.  The  number  of  design  variables  increased  from 
thirty-six  to  sixty-four  due  to  the  layered  quadrilateral  membrane  elements. 
Table  11  lists  material  properties,  minimum  sizes  constraints  and  fiber 
orientation.  The  NASTRAN  compatible  input  data  is  displayed  in  Figure  23. 

Results  of  Case  1  analysis  reside  in  Table  12  and  Figure  24. 

The  solution  procedure  path  followed  in  OPTFORCE  II  is  depicted  below: 

•  Initial  guess  for  design  variable  vector:  Minimum  size  constraint 
with  the  stress  ratio  option,  *  107.40  lbs. 

•  Program  entered  the  Linear  Programming  phase:  "Solution  Not 
Feasible";  Use  fully  stressed  design  (FSD)  as  minimum  weight 
solution.  wpgD  ■  110.73  lbs. 

•  Performed  check  on  the  Lagrange  multiplier  (X.)  calculations 
which  used  the  design  variable  vector  from  FSd.  Since  all 

X  >0  the  optimization  routine  terminated  with  the  FSD.  Note 
lthat  in  this  application  there  are  thirty-six  X's  associated 
with  the  design  variables,  thirty-six  X's  associated  with 
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TABLE  10  MATERIAL  PROPERTIES  A  CONSTRAINTS  -  CASE  I  SWEPT  WINGBGX 


(1)  Material  Properties:' 

Aluminum: 

E  -  10.6  x  106  psi 
p  »  .10  lbs/in* 
v  -  .25 

(2)  Minimum  sizes  (size  constraints) 

Quadrilateral  membranes  t  -  .10  in. 

*  2 

Rods  A  «*  .05  in. 

Shear  webs  t  •  ,05  in. 

w 

(3)  Allowable  stresses  (stress  constraints) 

o,  „  -  -30,000.0  psi,  o  -  30,000  psi 
lower  ’  r  upper 
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Figure  22  Optforce  II  Input  Date  -  Case  1  Swept  Wing box 


107 


TABLE  11  MATERIAL  PROPERTIES  &  CONSTRAINTS  -  CASE  2  SWEPT  VINGBOX 


(1)  Material  Properties: 

Aluainua:  E  -  10.5  x  10s  psi,  v  -  .30,  p  K  .10  lbs/in3 

Graphite/Epoxy:  Eu  -  18.5  x  10*  psi,  E22  *  1.6  x  10*  psi 
G  *  .65  x  10*  psi,  Vi2  -  .208,  V21  ■  .0203 
p  *  .055  lbs/lns 

(2)  Minimum  sizes  (size  constraints) 

Quadrilateral  membranes  t  -  .025  in. 

Rods  A*-  .05  in2 

Shear  webs  t  *  .050  in 

w 

(3)  Allowable  stresses  (stress  constraints) 

Aluminum  0.  -  -30000.0  psi,  O  *  30000.0  psi 

lower  r  upper  r 

Graphite/Epoxy  °iower  *  -110000.0  psi,  0upper  *  110000  psi 


(4)  Fiber  orientation 

Layer  No.  1  (top)  0  -  0° 

Layer  No.  2  0  -  45° 

Layer  No.  3  0  -  -45° 

Layer  No.  4  (bottom)  0-90 
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Figure  23  Opt force  II  Input  Data 
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TABLE  12  SWEPT  WINGBOX  RESULTS  -  CASE  1,  ALUMINUM  MATERIAL 


El. 

No. 

Element 

Type 

H 

El. 

No. 

Element 

Type 

H 

k* 

Quad 

.1000 

27109.0 

22 

Rod 

.0500 

198,9.0 

Quad 

.1000 

27101.0 

23 

Rod 

.0500 

13616.0 

Quad 
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Rod 
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Rod 
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29946.0 

29 

Rod 

.0500 

15907.0 

9 

Web 

.0500 

24010.0 

30 

Rod 

.0500 

6154.0 

10 

Web 

.1611 

30000.0 

31 

Rod 

.0500 

7886.0 

11 

Web 

.0587 

30000.0 

32 

Rod 

.0500 

4767.0 

12 

Web 

.0500 

6874.0 

33 

Rod 

.0500 

6235.0 

13 

Web 

.0500 

17029.0 

34 

Rod 

.0500 

3639.0 

14 

Web 

.0500 

7287.0 

35 

Rod 

.0500 

2792.0 

15 

Web 

.0500 

9057.0 

36 

Rod 

.0500 

2081.0 

16 

Web 

.1211 

30000.0 

W< 

17 

Web 

.0500 

17783.0 

Wa  -  110.73  lbs. 

18 

Web 

.0742 

30000.0 

Zio  »  12.63  in. 

19 

Web 

.0500 

1684.0 

CPU  -  24.08  sec 

.** 

20 

Web 

.0500 

27382.0 

I  -  1 

21 

Web 

.0500 

8031.0 

*A^  »  Design  variable  value,  rod  cross-sectional  area  (in2)  web  thickness 
(in),  quad  membrane  thickness  (in) 


**  IBM  3033 

i  ^ 

+  Stress  constraint  quantity:  ■  A -  1  4  0 

where  <Jj*  is  the  yield  stress  or  some  other  failure  stress  for  the  i— 
finite  element. 
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thirty-six  X's  associated  with  internal  stress  fields  and 
thirty-three  X's  associated  with  the  redundants • f or  a  total 
of  10S.  The  design  variable  and  stress  X's  sub  to  the  minimum 
weight  value  of  V  <*110.73  lbs.  This  Is  a  verification  of  the 
theoretical  development  given  in  Section  2.0. 


•  The  final  design  variable,  displacement,  element  stress  and 
reaction  vectors  are  calculated  and  presented  to  the  User  for 
review.  Minimum  weight  value  is  W  ■  110.73  lbs. 

A 

Examination  of  Table  12  and  Figure  24  shows  the  design  variable  vector 
and  element  stress  constraint  quantity  Y^/A^.  The  initial  and  final  weight 
of  the  structure  is  also  depicted  with  a  maximum  vertical  displacement  of 
12.63  Inches.  IBM  3033  computer  time  was  recorded  to  be  24.08  seconds. 

The  number  of  iterations,  1*1,  is  tabulated  to  indicate  the  Linear  Program¬ 
ming  phase  operation.  Figure  24  displays  the  design  variable  vector  in 
pictorial  form. 

The  solution  procedure  for  Case  2  was  as  follows: 

•  Initial  guess  vector:  Minimum  size  constraint  with  stress  ratio 
option,  W ^  •  60.05  lbs. 

•  Program  entered  the  Linear  Programming  phase:  "Solution  Mot 
Feasible";  Use  fully  stressed  design  as  minimum  weight  solution, 

W  -  83.32  lbs. 

IA 

•  Performed  check  on  Lagrange  Multiplier  (X  )  calculation.  Since 
certain  X's  associated  with  the  design  variables  were  negative 
the  full  Mewton-Raphson  routine  was  entered.  This  particular 
application  of  0PTF0RCE  II  required  fifty-four  X's  associated 
with  the  design  variable  vector,  fifty-four  X's  associated  with 
internal  stress  distribution  and  one-hundred  twenty-three  X's 
associated  with  the  redundants  for  a  total  of  231. 

•  Solution  procedure  terminated  after  four  iterations  in  the 
full  M-R  routine  due  to  an  excessive  numerical  increment 
associated  with  the  first  redundant  causing  structural  weight 
divergence.  User  nay  use  the  FSD  solution. 

Table  13  shows  pertinent  data  obtained  from  the  above  solution 
procedure.  Mote  the  drop  in  minimum  weight  through  use  of  the  graphite/ 
epoxy  wingbox  skins  as  compared  to  the  all  aluminum  wlngbox  structural 
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TABLE  13  SWEPT  WINGBOX  RESULTS  -  CASE  2,  GRAPHITE/ EPOXY  &  ALUMINUM  MATERIALS 


El. 

No. 

Element 

Type 

A  * 

Ai 

El. 

No. 

Element 

Type 

* 

Ai 

*i  + 

1 

Quad 

.0250 

41417.0 

31 

Rod 

.0500 

6560.0 

2 

Quad 

.0250 

48073.0 

32 

Rod 

.0500 

3308.0 

3 

Quad 

.0250 

41611.0 

33 

Rod 

.0500 

8279.0 

A 

Quad 

.0250 

41496.0 

34 

Rod 

.0500 

4512.0 

5 

Quad 

.0250 

11889.0 

35 

Rod 

.0500 

1219.0 

6 

Quad 

.0250 

18946.0 

36 

Rod 

.0500 

5853.0 

7 

Web 

.0655 

30000.0 

37 

Quad 

.0250 

62091.0 

8 

Web 

.0875 

30000.0 

38 

Quad 

.0250 

67220.0 

9 

Web 

.0520 

29968.0 

39 

Quad 

.0250 

51672.0 

10 

Web 

.1654 

30000.0 

40 

Quad 

.0250 

33913.0 

11 

Web 

.0500 

27775.0 

41 

Quad 

.0250 

26318.0 

12 

Web 

.0500 

3625.0 

42 

Quad 

.0250 

11784.0 

13 

Web 

.0500 

13716.0 

43 

Quad 

.0250 

3079.0 

14 

Web 

.0500 

5412.0 

44 

Quad 

.0250 

10121.0 

15 

Web 

.0500 

12279.0 

45 

Quad 

.0250 

4762.0 

16 

Web 

.1368 

30000.0 

46 

Quad 

.0250 

18134.0 

17 

Web 

.0500 

14070.0 

47 

Quad 

.0250 

23628.0 

18 

Web 

.0784 

30000.0 

48 

Quad 

.0250 

2260.0 

19 

Web 

.0500 

1374.0 

49 

Quad 

.0250 

6327.0 

20 

Web 

.0500 

21784.0 

50 

Quad 

.0250 

4139.0 

21 

Web 

.0500 

10887.0 

51 

Quad 

.0250 

14479.0 

22 

Rod 

.0500 

17595.0 

52 

Quad 

.0250 

10712.0 

23 

Rod 

.0500 

18751.0 

53 

Quad 

.0250 

10238.0 

24 

Rod 

.0500 

6953.0 

54 

Quad 

.0250 

10683.0 

TC 

Rod 

4J 

JUUUU  •  u 

26 

Rod 

.7203 

30000.0 

«i 

-  60.05  lbs. 

27 

Rod 

.0500 

12794.0 

wm 

-  83.32  lbs. 

28 

Rod 

.0500 

24769.0 

Z\  D 

-  13.94  in. 

29 

Rod 

.0500 

16438.0 

CPU 

*  249.54  sec** 

30 

Rod 

.0500 

10148.0 

I 

-  5 

*A^  •  Design  variable  value,  rod  **  IBM  3033 

cross-sectional  area  (in2), 

web  thickness  (in),  quad  +  Stress  constraint  quantity 

membrane  thickness  (in) 


weight.  Computer  time  was  registered  to  be  249.54  cpu  seconds  almost  a 
ten-fold  Increase  over  the  all-metallic  wing  solution  time.  This  is  a 
significant  amount  that  the  User  should  note  for  future  applications  of 
OPTFORCE  II.  The  increase  in  cpu  time  can  be  attributed  to  the  increased 
number  of  design  variables  and  Lagrange  multipliers.  One  hundred  and 
five  X's  were  needed  for  the  aluminum  wing  box  solution.  This  number 
increased  to  two  hundred  and  thirty-one  for  the  composite  wing  solution. 
The  number  of  iterations  1*5  is  tabulated  to  indicate  the  Linear  Program¬ 
ming  and  full  Newton-Raphson  phase  computations. 

It  is  of  interest  to  note  that  two  additional  analyses  were  con¬ 
ducted  to  show  the  effect  of  convergence  criteria  (OPTIM  input  card)  and 
minimum  size  constraint  imposed  on  the  graphite/epoxy  quadrilateral 
membrane  elements  (PQDMEM1  input  card) .  The  initial  convergence  criteria 
of  .010  was  changed  to  .0010.  Results  obtained  were  essentially  those 
shown  in  Table  13  except  that  the  full  N-R  iteration  routine  was 
terminated  due  to  an  excessive  delta  thickness  experienced  on  element 
number  18.  Thus,  in  this  particular  application  the  change  in  converg¬ 
ence  criteria  had  very  little  effect  upon  determining  the  minimum  weight 
design  with  the  exception  of  increasing  cpu  time  from  249.54  sec  to 
287.41  sec.  The  second  additional  analysis  completed  reduced  the  minimum 
thickness  constraint  of  the  membrane  elements  from  .0250  inch  to  .0100 
inches.  The  solution  procedure  again  followed  those  steps  initially 
shown.  The  minimum  weight  value  dropped  to  73.47  lbs.  as  would  be 
expected  and  the  deflection  of  the  wingbox  decreased  slightly  to  13.49 
inches.  Cpu  time  increased  about  45  seconds  to  a  value  of  290.06.  The 
distribution  of  the  design  variables  for  this  case  is  shown  in  Figure  25. 
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Y 


1 


Quad,  Membrane  Thickness 
Distribution 


tmin 


.025  in. 


*  all  layers  have  equal 
thickness 


(b)  Web  thickness  Distribution 


min 


.050  in. 


(.050) 

(.050) 


(c)  Rod  Area  Distribution 

A  .  -  .050  in? 

min 


(.050) 

(.050) 


Figure  25  Design  Variable  Distribution  -  Swept  Wingbox 

-  Case  2,  Graphite-Epoxy  &  Aluminum  Materials,  W  -73.47  lbs. 

n 
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4.0  CONCLUSIONS  &  RECOMMENDATIONS 


The  pilot  optimization  program,  OPTFORCE  I  (Ref.  4),  has  been 
succeaafully  extended  to  include  new  finite  elements,  optimization  con¬ 
straints  and  analysis  type  capabilities.  Additionally,  it  has  been  expanded 
into  a  general  purpose  type  optimization  code  featuring  NASTRAN  compatible 
input  data  formats  and  engineering  User  output  features.  As  a  result,  it 
is  concluded  that  the  technical  requirements  stated  in  Section  C  of  the 
subject  contract  F33615-80-C-3214  documents  have  been  met.  This  has 
resulted  in  a  new  optimization  code  labeled  OPTFORCE  II,  the  mathematical 
basis  of  which  has  been  amply  described  in  this  volume.  The  use  of  this 
code  has  also  been  illustrated  herein  and  its  input/output  features  and 
programming  aspects  are  given  in  Volume  Two  of  this  publication.  Ref.  11. 

Specific  tasks  which  came  to  a  successful  conclusion  were  the 
formulation  of  finite  element  matrices  bas-d  upon  the  force  method  approach. 
Membrane  triangle,  membrane  quadrilateral,  shear  panel  and  bar  (axial  force) 
elements  were  developed.  The  resultant  formulations  proved  to  be  accurate 
for  use  In  the  prediction  of  both  static  and  dynamic  behavior  of  structural 
components.  Their  use  in  optimization  analyses  also  proved  to  be  success¬ 
ful.  Formulation  of  optimization  constraint  equations  were  straight¬ 
forward  and  provided  for  the  first  time  the  opportunity  to  optimize 
structures  including  variable  stress,  multiple  displacement,  maximum  and 
minimum  size  and  dynamic  stiffness  (natural  frequency)  constraints  all 
within  the  context  of  the  force  method.  Multiple  load  conditions  are  also 
Included  in  the  formulation. 

The  force  method  formulation  of  a  rapid  re-analysis  technique 
concluded  in  a  highly  efficient  means  for  evaluating  the  effect  of  damage 
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on  the  stacic  and  dynamic  response  of  optimized  structures.  The  general 
theory  developed  in  Section  2.4  Included  the  use  of  the  aforementioned 
finite  elements,  however,  resources  only  permitted  the  coding  of  a  rapid 
re-analysis  program  using  bar  elements.  This  code  is  described  and 
illustrated  in  Appendix  B  of  Volume  II  and  is  successfully  demonstrated 
in  Section  2.4  of  this  volume.  Its  potential  is  proven  therein  by  the 
applications  shown  and  provides  a  means  for  determining  the  residual  strength 
and  response  of  damaged  structures.  The  method  is  a  direct  one  and  elimin¬ 
ates  the  Iterative  nature  of  other  methods  presently  available  in  the 
literature. 

Executions  of  OPTFORCE  II  were  ample  enough  to  conclude  that  the 
initial  capabilities  of  the  force  method  code  OPTFORCE  I  could  be  success¬ 
fully  extended  to  optimize  three-dimensional  truss-like  and  aerospace  type 
structures.  The  applications  of  OPTFORCE  II  further  proved  that  the  force 
method  yields  more  accurate  optimization  solutions  than  the  displacement 
based  OPTIM  III  reference  computer  code.  The  efficiency  studies  conducted 
and  profusely  described  in  Section  3.1  showed  this  fact  to  be  true.  These 
studies  further  demonstrated  that  computer  time  (epu  seconds)  is  a  less 
relevant  measure  of  computer  code  "efficiency"  than  originally  thought. 

It  is  concluded  from  these  studies  that  the  OPTFORCE  II  code  is  preferred 
over  the  reference  displacement  based  OPTIM  III  code. 

No  difficulties  were  experienced  with  the  use  of  OPTFORCE  II  when¬ 
ever  structures  were  optimized  subject  to  minimum  size  and  stress  constraints. 
Minimum  weight  solutions  including  displacement  and/or  dynamic  constraints 
in  addition  to  size  and  stress  constraints  proved  troublesome.  No  satis¬ 
factory  solutions  could  be  found  for  the  numerical  difficulties  encountered. 
Obviously  additional  resources  are  needed  to  obviate  this  problem.  It  is 


concluded  at  this  tine  that  the  difficulty  lies  not  in  the  force  method 
foraulation  of  those  type  of  constraints  but  in  the  numerical  application 
of  attendant  governing  equations.  It  is  recommended  that  further  studies 
of  this  perplexity  be  pursued. 

Finally,  it  is  concluded  that  the  research  conducted  and  reported 
herein  has  shown  that  the  force  method  of  structural  optimization  has  pro¬ 
vided  the  solution  to  the  technical 'problems  discussed  in  the  Introduction 
to  this  report.  The  force  method  approach  gives  the  means  to  express 
constraints,  in  particular  the  stress  constraint,  in  such  a  manner  as  to 
provide  the  needed  mathematical  expressions  for  use  in  the  Lagrangan 
formulation  of  the  weight  optimization  solution  procedure.  This  enhances 
the  definition  of  the  optimality  criteria  and  attendant  first  and  second 
derivatives  used  in  the  nonlinear  solution  procedure.  Thus ,  it  is  further 
recommended  that  the  force  method  optimization  procedure  resident  in 
OPTFQRCE  II  be  seriously  considered  as  the  future  general  purpose  struc¬ 
tural  optimization  code. 
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